Actividad 3: Regresión Lineal Múltiple
Víctor Escandón
2023-09-29
Ejercicio 1
Sesgo por variables omitidas
library(pacman)
p_load(data.table, fixest, magrittr, sandwich, mvtnorm, MASS)
set.seed(123)# Load data
library(wooldridge)
data("bwght")
bwght$cigs2 <-bwght$cigs*bwght$cigs #adding a column with squares
# Run regressions
est_0 =feols(bwght ~ cigs, data = bwght)
est_1 = feols(bwght ~ cigs + motheduc, data = bwght)
est_2= feols(bwght ~ cigs + motheduc+faminc, data = bwght)
est_3 = feols(bwght ~ cigs + cigs2 + motheduc + faminc, data = bwght)#Nicer tables
## Standard errors (using fixest package)
etable(est_0, est_1, est_2, est_3, se = "standard") # Under MLR5## est_0 est_1 est_2
## Dependent Var.: bwght bwght bwght
##
## Constant 119.8*** (0.5723) 115.4*** (3.107) 116.8*** (3.138)
## cigs -0.5138*** (0.0905) -0.4862*** (0.0926) -0.4633*** (0.0927)
## motheduc 0.3308 (0.2328) 0.0143 (0.2580)
## faminc 0.0915** (0.0325)
## cigs2
## _______________ ___________________ ___________________ ___________________
## S.E. type IID IID IID
## Observations 1,388 1,387 1,387
## R2 0.02273 0.02420 0.02977
## Adj. R2 0.02202 0.02279 0.02767
##
## est_3
## Dependent Var.: bwght
##
## Constant 117.3*** (3.150)
## cigs -0.7787*** (0.2099)
## motheduc -0.0074 (0.2581)
## faminc 0.0898** (0.0324)
## cigs2 0.0123. (0.0073)
## _______________ ___________________
## S.E. type IID
## Observations 1,387
## R2 0.03174
## Adj. R2 0.02894
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1Usando est_2 prediga el valor en onzas del peso de un niño cuya madre fumó 2 cigarros por día durante el embarazo, con todo lo demás en su promedio.
En promedio, si la madre fumó dos cigarros, entonces el peso del bebe es afectado dos veces el estimador de cigarros sobre preso de bebe, ceteris paribus. Es decir, \(2*-0.4633=-0.9266\)
¿Cuál es el signo de la correlación entre educación de la madre y número de cigarros fumados al día durante el embarazo? ¿Es lo que esperaba? Explique. (Vea Tabla 3.2. en Wooldridge ed. 7)
bwght_cleaned <- na.omit(bwght[, c("cigs", "motheduc")])
correlation <- cor(bwght_cleaned$cigs, bwght_cleaned$motheduc)
print(correlation)## [1] -0.2138651La correlación es negativa. Tiene sentido, ya que el estimador de cigarros sobre el peso del bebé es menor cuando no tomamos en cuenta la educación de la madre. Cuando se toma en cuenta, dicho estimador es un poco más alto. Esto debido a que el estimador estaba siendo infraestimado por omitir la variable de la educación de la madre.
- ¿Cuál es el efecto de fumar un cigarro extra en el embarazo en el peso en onzas considerando una relación no lineal? ¿Por qué utilizaría una relación no lineal? Explique.
El efecto está dado por \[\frac{\partial\hat{bwght}}{\partial\hat{cigs}}= \beta_1+2\beta_2\hat{cigs}\] Es decir, en promedio, el aumento del consumo de un cigarro durante el embarazo es de \(-0.7787+2(0.0123)=-0.7541\).
- Interprete R2 y R2-adj en est_3. Formalmente ¿Por qué R2-adj es menor?
\(R^2\) expresa que las 4 variables utilizadas en la reg3 explicaran 0.03174 la varianza del peso del bebè, no obtante, este valor es demasiado alto, ya que al insertar más regresores, R2 por definición aumenta. \(adjR^2\) corrige esta alteración. En realidad, las variables utilizadas en la regresión realmente explican 0.02894 la variación de la variable dependiente.
Ejercicio 2
## Standard errors (using fixest package)
#Table 1:
etable(est_0, est_1, est_2, est_3, se = "standard") # Under MLR5## est_0 est_1 est_2
## Dependent Var.: bwght bwght bwght
##
## Constant 119.8*** (0.5723) 115.4*** (3.107) 116.8*** (3.138)
## cigs -0.5138*** (0.0905) -0.4862*** (0.0926) -0.4633*** (0.0927)
## motheduc 0.3308 (0.2328) 0.0143 (0.2580)
## faminc 0.0915** (0.0325)
## cigs2
## _______________ ___________________ ___________________ ___________________
## S.E. type IID IID IID
## Observations 1,388 1,387 1,387
## R2 0.02273 0.02420 0.02977
## Adj. R2 0.02202 0.02279 0.02767
##
## est_3
## Dependent Var.: bwght
##
## Constant 117.3*** (3.150)
## cigs -0.7787*** (0.2099)
## motheduc -0.0074 (0.2581)
## faminc 0.0898** (0.0324)
## cigs2 0.0123. (0.0073)
## _______________ ___________________
## S.E. type IID
## Observations 1,387
## R2 0.03174
## Adj. R2 0.02894
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1#Table 2:
etable(est_0, est_1, est_2, est_3, se = "hetero") # Fixes Heteroskedasticity (HC1 type) ## est_0 est_1 est_2
## Dependent Var.: bwght bwght bwght
##
## Constant 119.8*** (0.5745) 115.4*** (2.928) 116.8*** (2.945)
## cigs -0.5138*** (0.0877) -0.4862*** (0.0890) -0.4633*** (0.0893)
## motheduc 0.3308 (0.2202) 0.0143 (0.2408)
## faminc 0.0915** (0.0312)
## cigs2
## _______________ ___________________ ___________________ ___________________
## S.E. type Heteroskedast.-rob. Heteroskedast.-rob. Heteroskedast.-rob.
## Observations 1,388 1,387 1,387
## R2 0.02273 0.02420 0.02977
## Adj. R2 0.02202 0.02279 0.02767
##
## est_3
## Dependent Var.: bwght
##
## Constant 117.3*** (2.962)
## cigs -0.7787*** (0.2018)
## motheduc -0.0074 (0.2410)
## faminc 0.0898** (0.0312)
## cigs2 0.0123. (0.0070)
## _______________ ___________________
## S.E. type Heteroskedast.-rob.
## Observations 1,387
## R2 0.03174
## Adj. R2 0.02894
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1- Demuestre matemáticamente por qué el SE en la segunda tabla es distinto al SE la Tabla 1, al menos para el estimador asociado a cigs. (i.e derive la varianza del error homoscedastica y heteroscedástica).
Dado que para para el caso homosedástico \[{Var}(\hat{\beta}_j)= \frac{\sigma^2} {Var(x)}\] y para el caso heterosedástico \[{Var}(\hat{\beta}j)= \frac{\sigma^2\hat{u^2_{j}}}{Var(x)}\] es notorio que la varianza está alterada en razón de una situación heterosedástica u homosedástica. Por tanto, el error estándar es diferente para ambos casos.
- ¿Qué supuesto asumimos como vulnerado en la estimación de Tabla 2?
El supuesto MLS.5. La varianza de los errores debe ser constante. \[Var(u_i|\textbf{x}_i)=\sigma^2\]
- Dado:
summary(lm(bwght$cigs ~ bwght$motheduc))y considerando est_0 (en ejercicio 1) ¿aumenta la varianza del estimador cigs al incluirse mothereduc en est_1? Discuta lo anterior utilizando las fórmulas de la varianza de beta para un modelo simple y para uno múltiple.
summary(lm(bwght$cigs ~ bwght$motheduc)) ##
## Call:
## lm(formula = bwght$cigs ~ bwght$motheduc)
##
## Residuals:
## Min 1Q Median 3Q Max
## -7.968 -2.592 -2.054 -0.441 47.408
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 9.04312 0.86783 10.420 < 2e-16 ***
## bwght$motheduc -0.53761 0.06598 -8.148 8.21e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 5.838 on 1385 degrees of freedom
## (1 observation deleted due to missingness)
## Multiple R-squared: 0.04574, Adjusted R-squared: 0.04505
## F-statistic: 66.38 on 1 and 1385 DF, p-value: 8.212e-16Dado que las formulas de la varianza de los estimadores (bajo homosedasticidad) son:
Regresión simple \[Var(\hat{\beta_j})=\frac{\sigma^2}{SST_x^2}\]
Regresión múltiple: \[Var(\hat{\beta_j})=\frac{\sigma^2}{SST_x^2(1-R_j^2)}\]
Notamos que, agregar una variable nueva no afecta la varianza cuando se analiza una regresión simple. No obstante, lo correcto es analizar la varianza en el caso donde ahora se realiza una regresión múltiple. Notamos que, ahora, la nueva variable afecta la varianza del estimador de cigs. Específicamente, \(R_j^2\) proviene de una regresión de cigs sobre motheduc. Por tanto, note que al agrergar un regresor implica agregar dicho nuevo término que hará que la variación cigs aumente.
- ¿Cuál es el VIF Tras agregar motheduc? Discuta, considerando el trade-off sesgo varianza si debemos incluir motheduc.
Para calcular VIF, primero es necesario obtener \(R^2\) de correr cigs sobre motheduc.
regcigsmotheduc <- lm(cigs ~ motheduc, data = bwght, na.action = na.exclude)
r2cigsmotheduc<-summary(regcigsmotheduc)[8]
print(r2cigsmotheduc)## $r.squared
## [1] 0.04573828Así, VIF de cigs tras agregar motheduc es:
VIFcigs<-1/(1- 0.04573828)
print(VIFcigs)## [1] 1.047931La discusión en la siguiente: agregar un regresor implica aumento a la varianza del estimador de, en este caso, cigs, aunque también implica menos sesgo. Quiere decir que aumenta la probabilidad de no obtener al valor poblacional, pero se estará cerca de él. Por otro lado, no agregar dicho regresor implica menos varianza, pero más sesgo. Esto significa que el valor de la estimación será poco invariable, pero nunca cerca del valor poblacional. Por ello, siempre es preferible tener menos sesgo y más varianza. Además, la varianza puede ser resuelta aumentando el número de observaciones.
- ¿A qué concepto nos referimos en el inciso c?
El problema de colinearidad
Ejercicio 3
The data BEAUTY (Wooldridge Library in R) contains data on wages and an index of people’s beauty.
- Explore las variables con ?beauty y agregue una columna a los datos
que incluya una aproximación de edad para cada i:
beauty$age <- beauty$educ+beauty$exper+6
beauty$age <- beauty$educ+beauty$exper+6- Estime una regresión simple: de log(wage) con looks como explicativa con errores estándar robustos a heteroscedasticidad (HC1): pegue el resultado abajo e interprete B1
reg3.b<-feols(lwage~looks,data=beauty)
etable(reg3.b,se="white")## reg3.b
## Dependent Var.: lwage
##
## Constant 1.496*** (0.0818)
## looks 0.0510* (0.0252)
## _______________ _________________
## S.E. type Heteroskeda.-rob.
## Observations 1,260
## R2 0.00345
## Adj. R2 0.00265
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1En promedio, aumentar en una unidad el look, aumenta 5.1% el salario, ceteris paribus.
- Agregue como controles a la regresión en b. la experiencia, la experiencia al cuadrado, la educación en años y la edad ¿Qué pasa con su regresión y por qué? ¿Qué supuesto se puede estar violando?
reg3.c<-feols(lwage~looks+exper+expersq+educ+age,data=beauty)## The variable 'age' has been removed because of collinearity (see $collin.var).Se elimió la variable age, esto porque viela es supuesto de no colinearidad. R automáticamente elimina una de las variables; en este caso, age.
- Solucione el problema en c. y vuelva a estimar su regresión. Pegue el resultado abajo e interprete B1.
etable(reg3.c)## reg3.c
## Dependent Var.: lwage
##
## Constant 0.0463 (0.1053)
## looks 0.0621** (0.0219)
## exper 0.0487*** (0.0048)
## expersq -0.0007*** (0.0001)
## educ 0.0684*** (0.0058)
## _______________ ___________________
## S.E. type IID
## Observations 1,260
## R2 0.23225
## Adj. R2 0.22981
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1R elimina automáticamente la variable age, es decir, ya soluciona el problema. Así que, en promedio, al aumentar en una unidad el look, el salario aumenta 6.2%, ceteris paribus.
Ejercicio 4
Which of the following can cause OLS estimators to be biased? Why? a. Heteroskedasticity. b. Omitting an important variable. c. A sample correlation coefficient of .95 between two independent variables both included in the model.
La respuesta es b, ya que se obtiene el sesgo por variable omitida
Ejercicio 5
- Which of the following can cause the variance of OLS estimators to be biased? Why?
- Heteroskedasticity.
- Omitting an important variable.
- A sample correlation coefficient of .95 between two independent variables both included in the model.
La respuesta es a. Una variación no constante de los errores implica un sesgo en la varianza del estimador.
Ejercicio 6
- Confirm the partialling out interpretation of the OLS estimates by explicitly doing the partialling out for Example 3.2 (see Wooldridge 7th Edition Section 3.2f). This first requires regressing educ on exper and tenure and saving the residuals, r1. Then, regress log(wage) on r1. Compare the coefficient on r1 with the coefficient on educ in the regression of log(wage) on educ, exper, and tenure.
regPO<-feols(educ~exper+tenure,data=wage1)
r1<-regPO$residuals
regPO.1<-feols(lwage~r1,data=wage1)- Pegue abajo el resultado de las dos regresiones de interés y explique intuitivamente por qué se observa que el coeficiente de r1 es igual al de educ cuando controla por exper y tenure.
summary(regPO.1)## OLS estimation, Dep. Var.: lwage
## Observations: 526
## Standard-errors: IID
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 1.623268 0.020664 78.5553 < 2.2e-16 ***
## 0.092029 0.007880 11.6794 < 2.2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## RMSE: 0.473021 Adj. R2: 0.205037reg6<-feols(lwage~educ+exper+tenure,data=wage1)
summary(reg6)## OLS estimation, Dep. Var.: lwage
## Observations: 526
## Standard-errors: IID
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 0.284360 0.104190 2.72923 6.5625e-03 **
## educ 0.092029 0.007330 12.55525 < 2.2e-16 ***
## exper 0.004121 0.001723 2.39144 1.7136e-02 *
## tenure 0.022067 0.003094 7.13307 3.2944e-12 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## RMSE: 0.439183 Adj. R2: 0.312082Para ambas el coeficiente es 0.092029. Esto ya que se aplicó al teorema de Frisch-Waugh. Esto ya que la primera regresión, se corrió una regresión para mandar los efectos de educ sobre wage a los residuos. Esto hace que correr una regresión sobre dichos residuos me arroja lo mismo que correr una regresión tomando en cuenta la variable de educ.
Ejercicio 7
# set seed
set.seed(123)
#Example before computation:
#we define a matrix X with two vectors (regressors x1 and x2) from random normal values with mean [5,10]
#and a matrix "sigma" of variances [3,6], and covariance 2
mu <- c(5,10)
sigma <- matrix(c(3,2,2,6), 2, 2)
Xej7 <- rmvnorm(1000, mean = mu, sigma = sigma)
?rmvnorm
head(Xej7)## [,1] [,2]
## [1,] 3.955929 9.171603
## [2,] 7.623002 10.937024
## [3,] 6.059552 14.178891
## [4,] 5.142173 7.191611
## [5,] 3.639927 8.592305
## [6,] 7.209863 11.466354summary(Xej7)## V1 V2
## Min. : 0.06447 Min. : 3.454
## 1st Qu.: 3.92522 1st Qu.: 8.403
## Median : 5.16410 Median :10.059
## Mean : 5.07445 Mean :10.066
## 3rd Qu.: 6.11615 3rd Qu.:11.795
## Max. :11.04032 Max. :19.072var(Xej7)## [,1] [,2]
## [1,] 2.768283 1.727895
## [2,] 1.727895 6.114899#plot vectors 1 and 2 in our matrix X (with covariance 2!)
plot(Xej7[,1], Xej7[,2])
##end of example.# Set number of observations
n <- 500
#Another example of what is going to happen... Here we create a matrix X
#with two vectors (regressors x1 and x2) with median 50 and 100 each, and cov(X_1,X_2) = 0.25
Xej7.2 <-rmvnorm(n, c(50, 100), sigma = cbind(c(10, 2.5), c(2.5, 10)))
summary (Xej7.2)## V1 V2
## Min. :41.58 Min. : 92.43
## 1st Qu.:47.90 1st Qu.: 97.96
## Median :50.00 Median : 99.71
## Mean :49.98 Mean : 99.88
## 3rd Qu.:51.98 3rd Qu.:102.04
## Max. :59.31 Max. :110.06plot(Xej7.2[,1], Xej7.2[,2]) 
#Now the loop to do this 10000 times, fasten your seat belts:
# Empty vectors of coefficients to be filled
coefs1 <- cbind("beta_hat_1" = numeric(10000), "beta_hat_2" = numeric(10000))
coefs2 <- coefs1
###From here###
# Loop sampling and estimation (be patient it takes a while)
for (i in 1:10000) {
# for median 50 and 100 and cov(X_1,X_2) = 0.25
Xej7.2 <- rmvnorm(n, c(50, 100), sigma = cbind(c(10, 2.5), c(2.5, 10)))
uej7.2 <- rnorm(n, sd = 5)
#We define PRF
Yej7.2 <- 5 + 2.5 * Xej7.2[, 1] + 3.5 * Xej7.2[, 2] + uej7.2
#we compute SRF 10000 times using our n values in matrix X
coefs1[i, ] <- lm(Yej7.2 ~ Xej7.2[, 1] + Xej7.2[, 2])$coefficients[-1]
# for cov(X_1,X_2) = 0.85
Xej7.3 <- rmvnorm(n, c(50, 100), sigma = cbind(c(10, 8.5), c(8.5, 10)))
Yej7.3 <- 5 + 2.5 * Xej7.3[, 1] + 3.5 * Xej7.3[, 2] + uej7.2
coefs2[i, ] <- lm(Yej7.3 ~ Xej7.3[, 1] + Xej7.3[, 2])$coefficients[-1]
}
## Este código genera, con el uso de un loop, la linea poblacional que relaciona todos los valores de X con Y. Pero hace dos funciones de Y, una con valores de X con covarianza 2.5 y otra con cov 8.5 Después corre una regresión lineal que llena los valores del vector de beta.
#Histogram of coefficients x1 and x2 for the two 10000s loops of regressions.
hist(coefs1)
hist(coefs2)
## Generan histogramas de los coeficientes.
# Obtain variance estimates
diag(var(coefs1))## beta_hat_1 beta_hat_2
## 0.005266405 0.005323088diag(var(coefs2))## beta_hat_1 beta_hat_2
## 0.01817545 0.01825663# Obtiene varianzas de los coeficientes- A partir de la marca ###From here###. Explique cada línea pertinente del código
Vease el código
- Escriba la matriz de varianza – covarianza de los estimadores 1 y 2.
\[Var(\beta|X)=\sigma^2(\sum_{i=1}^{N}x_ix_i´)^{-1}\]
¿Por qué los histogramas de coefs1 y coefs2 son distintos?
¿Por qué las varianzas de coefs1 y coefs2 son distintas?
Porque las covarianzas en las X, para cada caso, tienen covarianza diferente.
- ¿qué parte del código genera estas diferencias?
Xej7.2 <- rmvnorm(n, c(50, 100), sigma = cbind(c(10, 2.5), c(2.5, 10)))
y
Xej7.3 <- rmvnorm(n, c(50, 100), sigma = cbind(c(10, 8.5), c(8.5, 10)))
- ¿Qué concepto/supuesto central estamos simulando en este ejercicio?
Multicolinearidad
- Ahora cambie n a 1000, ¿qué sucede con los histogramas y las varianzas? Explique intuitivamente, pero con el uso de alguna fórmula ¿por qué se da este cambio?
Dado que \[Var(\beta|X)=\sigma^2(\sum_{i=1}^{N}x_ix_i´)^{-1}\] y \(SST_x^2\) aumenta cuando n aumenta, entonces el denominador se hace más grande y, por tanto, la \(Var(\beta|X)\) disminuye. Esto hace que los histogramas se hagan más estrechos.
# Set number of observations
n.0 <- 1000
#Another example of what is going to happen... Here we create a matrix X
#with two vectors (regressors x1 and x2) with median 50 and 100 each, and cov(X_1,X_2) = 0.25
Xej7.2.0 <-rmvnorm(n.0, c(50, 100), sigma = cbind(c(10, 2.5), c(2.5, 10)))
summary (Xej7.2.0)## V1 V2
## Min. :39.32 Min. : 90.92
## 1st Qu.:47.94 1st Qu.: 97.94
## Median :49.97 Median :100.11
## Mean :50.03 Mean :100.08
## 3rd Qu.:52.15 3rd Qu.:102.11
## Max. :59.75 Max. :109.41plot(Xej7.2.0[,1], Xej7.2.0[,2]) 
#Now the loop to do this 10000 times, fasten your seat belts:
# Empty vectors of coefficients to be filled
coefs1.0 <- cbind("beta_hat_1" = numeric(10000), "beta_hat_2" = numeric(10000))
coefs2.0 <- coefs1.0
###From here###
# Loop sampling and estimation (be patient it takes a while)
for (i in 1:10000) {
# for median 50 and 100 and cov(X_1,X_2) = 0.25
Xej7.2.0 <- rmvnorm(n.0, c(50, 100), sigma = cbind(c(10, 2.5), c(2.5, 10)))
uej7.2.0 <- rnorm(n.0, sd = 5)
#We define PRF
Yej7.2.0 <- 5 + 2.5 * Xej7.2.0[, 1] + 3.5 * Xej7.2.0[, 2] + uej7.2.0
#we compute SRF 10000 times using our n values in matrix X
coefs1.0[i, ] <- lm(Yej7.2.0 ~ Xej7.2.0[, 1] + Xej7.2.0[, 2])$coefficients[-1]
# for cov(X_1,X_2) = 0.85
Xej7.3.0 <- rmvnorm(n.0, c(50, 100), sigma = cbind(c(10, 8.5), c(8.5, 10)))
Yej7.3.0 <- 5 + 2.5 * Xej7.3.0[, 1] + 3.5 * Xej7.3.0[, 2] + uej7.2.0
coefs2.0[i, ] <- lm(Yej7.3.0 ~ Xej7.3.0[, 1] + Xej7.3.0[, 2])$coefficients[-1]
}
## Este código genera, con el uso de un loop, la linea poblacional que relaciona todos los valores de X con Y. Pero hace dos funciones de Y, una con valores de X con covarianza 2.5 y otra con cov 8.5 Después corre una regresión lineal que llena los valores del vector de beta.
#Histogram of coefficients x1 and x2 for the two 10000s loops of regressions.
hist(coefs1.0)
hist(coefs2.0)
## Generan histogramas de los coeficientes.
# Obtain variance estimates
diag(var(coefs1.0))## beta_hat_1 beta_hat_2
## 0.002647207 0.002695637diag(var(coefs2.0))## beta_hat_1 beta_hat_2
## 0.009097768 0.009120378# Obtiene varianzas de los coeficientes- Aún bajo el supuesto MRL4 ¿Se puede obtener un estimador inconsistente en una muestra aleatoria de individuos? De un ejemplo con el uso de los histogramas.
Violando alguno de los otro 3 supuestos, se obtiene un paraámetro insesgado y, por tanto, inconsistente.