Python - Regresión Lineal en R

Regresión simple

modelo = lm(formula, data)
y_barra = mean(modelo$variabley)
x_barra = mean(modelo$variablex)
cov_xy = cov(modelo$variablex, modelo$y)
var_x = var(modelo$variablex)

hb1 = cov_xy/var_x
hb0 = y_barra-hb1*x_barra

Gráfica MCO

with(ceosal1,
plot(roe, salary,
ylim = c(0, 7500),
col = "dark gray",
main = "Salario y desempeño del CEO",
xlab = "Return on Equity (ROE)",
ylab = "Salario (miles de $USD)"))
# Grafica recta mco
abline(a = hb0, b = hb1,
col = "red", lwd = 1.5)

# Legend
legend(x = "topleft", box.col = "black",
lty = 1, col = "red",
legend = "Recta de MCO")

# Valores observados
with(wage1,
plot(educ, wage,
col = "black",
pch = 19, cex = 0.65,
xlim = c(0, 19),
ylim = c(0, 25),
main = "Educación y Salarios",
xlab = "Años de estudio",
ylab = "Salario (por hora)"))
legend(x = "topleft", box.col = "black",
pch = 19,
legend = c("Observaciones"))

# Valores ajustados y residuales
with(wage1,
plot(educ, hat_y,
col = "blue",
pch = 20, cex = 0.75,
xlim = c(0, 19),
ylim = c(0, 25),
main = "Valores ajustados y Residuales",
xlab = "Observaciones",
ylab = "Estadísticos"))
with(wage1,
points(educ, hat_u,
col = "cyan",
pch = 1, cex = 0.75))
legend(x = "topleft", box.col = "black",
pch = c(20, 1),
col = c("blue", "cyan"),
legend = c("Ajustados", "Residuales"))

Sesgo de variable omitida

# Estimar el modelo multiple
modelo_desempeno_academico <- lm(formula = colGPA ~ hsGPA + ACT,
data = gpa1)

# Estimar el modelo simple
modelo_desempeno_simple <- lm(formula = colGPA ~ hsGPA,
data = gpa1)

# Presentacion de los resultados
stargazer(modelo_desempeno_academico,
modelo_desempeno_simple,
digits = 3,
type = "text", style = "qje",
report = c("vc*s"),
column.labels = c("Múltiple", "Simple"),
dep.var.labels = c("Promedio Univ"),
covariate.labels = c("Promedio Prepa",
"Logros",
"Constante"),
keep.stat = c("n", "rsq"))
# Estimar el modelo simple del sesgo
modelo_sesgo <- lm(formula = ACT ~ hsGPA,
data = gpa1)

# Presentacion de los resultados
stargazer(modelo_sesgo, digits = 3,
type = "text", style = "qje",
report = c("vc*s"),
dep.var.labels = c("Logros"),
covariate.labels = c("Promedio Prepa",
"Constante"),
keep.stat = c("n", "rsq"))

Con esto, podemos calcular que el sesgo es:
β̂2δ1~=0.009×3.074=0.029
Más aún, podemos verificar la igualdad
β1~=β̂1+β̂2δ1

# Coeficiente sesgado estimado
round(modelo_desempeno_simple$coefficient[2], 3)

# Coeficiente sesgado teorico
round(modelo_desempeno_academico$coefficient[2] +
modelo_desempeno_academico$coefficient[3] * modelo_sesgo$coefficient[2], 3)

Problema de Multicolinealidad

# Estimacion del modelo
modelo_salarios <- lm(wage ~ educ + exper + tenure,
data = base)

# Presentacion de resultados
stargazer(modelo_salarios, digits = 3,
type = "text", style = "qje",
report = c("vc*s"),
dep.var.labels = c("Salario (por hora)"),
covariate.labels = c("Educación",
"Experiencia",
"Antigüedad",
"Constante"),
keep.stat = c("n", "rsq"))
# Calculo de los VIF
vifs <- vif(modelo_salarios)
vifs

Estos resultados indican que:

La varianza del estimador asociado a educ es 1.113 veces mayor respecto a su “mínimo posible”.
La varianza del estimador asociado a exper es 1.478 veces mayor respecto a su “mínimo posible”.
La varianza del estimador asociado a tenure es 1.349 veces mayor respecto a su “mínimo posible”.

Un VIF menor a cinco es considera como bueno.
Un VIF mayor a cinco, pero menor a diez es considerado no ideal, pero tampoco necesariamente un problema.
Un VIF mayor a 10 es considerado un problema y sugiere analizar la relevancia de incorporar la variable en el modelo.

Prueba t

#Código para prueba t cola derecha
prueba_t_rt <- function(modelo, variable, valor = 0, alpha = 0.05){

beta <- summary(modelo)$coefficients[variable,1]
se <- summary(modelo)$coefficients[variable,2]
t <- (beta-valor)/se


n <- nobs(modelo)
k <- nrow(summary(modelo)$coefficients)-1
df <- n-k-1

valcrit <- qt((1 - alpha), df)
p_val_der <- pt(-t,df)

if (t>valcrit){
cat("Se realizó una prueba de COLA DERECHA sobre el coeficiente de:",variable,
"La significancia de la prueba es del:",100*alpha,"%",
"La hipótesis nula a probar es H0:",variable,"tiene un efecto de valor:",valor,
"El estadístico t de la prueba es:",t,"y el valor crítico de la prueba es:",valcrit,
"El valor-p de la prueba es:",p_val_der,
"Por esto, se concluye que H0 SE RECHAZA")
}else{
cat("Se realizó una prueba de COLA DERECHA sobre el coeficiente de:",variable,
"La significancia de la prueba es del:",100*alpha,"%",
"La hipótesis nula a probar es H0:",variable,"tiene un efecto de valor:",valor,
"El estadístico t de la prueba es:",t,"y el valor crítico de la prueba es:",valcrit,
"El valor-p de la prueba es:",p_val_der,
"Por esto, se concluye que H0 NO SE RECHAZA")
}

}

#Código para prueba t cola izquierda
prueba_t_lt <- function(modelo, variable, valor = 0, alpha = 0.05){

beta <- summary(modelo)$coefficients[variable,1]
se <- summary(modelo)$coefficients[variable,2]
t <- (beta-valor)/se

n <- nobs(modelo)
k <- nrow(summary(modelo)$coefficients)-1
df <- n-k-1

valcrit <- qt(1 - alpha, df)
p_val_izq <- pt(t,df)

if (t < -valcrit){
cat("Se realizó una prueba de COLA IZQUIERDA sobre el coeficiente de:",variable,
"La significancia de la prueba es del:",100*alpha,"%",
"La hipótesis nula a probar es H0:",variable,"tiene un efecto de valor:",valor,
"El estadístico t de la prueba es:",t,"y el valor crítico de la prueba es:",valcrit,
"El valor-p de la prueba es:",p_val_izq,
"Por esto, se concluye que H0 SE RECHAZA")
}else{
cat("Se realizó una prueba de COLA IZQUIERDA sobre el coeficiente de:",variable,
"La significancia de la prueba es del:",100*alpha,"%",
"La hipótesis nula a probar es H0:",variable,"tiene un efecto de valor:",valor,
"El estadístico t de la prueba es:",t,"y el valor crítico de la prueba es:",valcrit,
"El valor-p de la prueba es:",p_val_izq,
"Por esto, se concluye que H0 NO SE RECHAZA")
}

}

#Código para prueba t dos colas
prueba_t <- function(modelo, variable, k, valor = 0, alpha = 0.05){

beta <- summary(modelo)$coefficients[variable,1]
se <- summary(modelo)$coefficients[variable,2]
t <- (beta-valor)/se

n <- nobs(modelo)
k <- nrow(summary(modelo)$coefficients)-1
df <- n-k-1

valcrit <- qt(1 - (alpha/2), df)
p_val <- 2*pt(-abs(t),df)

if (abs(t)>valcrit){
cat("Se realizó una prueba de DOS COLAS sobre el coeficiente de:",variable,
"La significancia de la prueba es del:",100*alpha,"%",
"La hipótesis nula a probar es H0:",variable,"tiene un efecto de valor:",valor,
"El estadístico t de la prueba es:",t,"y el valor crítico de la prueba es:",valcrit,
"El valor-p de la prueba es:",p_val,
"Por esto, se concluye que H0 SE RECHAZA")
}else{
cat("Se realizó una prueba de DOS COLAS sobre el coeficiente de:",variable,
"La significancia de la prueba es del:",100*alpha,"%",
"La hipótesis nula a probar es H0:",variable,"tiene un efecto de valor:",valor,
"El estadístico t de la prueba es:",t,"y el valor crítico de la prueba es:",valcrit,
"El valor-p de la prueba es:",p_val,
"Por esto, se concluye que H0 NO SE RECHAZA")
}

}

Prueba F

El modelo restringido es sin las variables que no se consideran importantes

# Estimacion del modelo restringido
modelo_baseball_res <- lm(formula = salary ~ years + gamesyr,
data = base3)

# Presentacion de los resultados
stargazer(modelo_baseball_nores,
modelo_baseball_res, digits = 3,
type = "text", style = "aer",
dep.var.labels = c("Salario"),
covariate.labels = c("Antigüedad", "Juegos/año",
"Promedio Bateo", "HR/año",
"RBIS/año", "Constante"),
keep.stat = c("n", "rsq", "f"))

# Elementos necesarios para el calculo del estadistico F
ssr_nores <- sum(modelo_baseball_nores$residuals^2)
ssr_res <- sum(modelo_baseball_res$residuals^2)
n <- nobs(modelo_baseball_nores)
k <- 5
q <- 3

# Calculo el estadistico F
f_stat <- ((ssr_res - ssr_nores) / ssr_nores) * ((n - k -1) / q)
f_stat

# Elementos necesarios para el calculo del estadistico F
r2_nores <- summary(modelo_baseball_nores)$r.squared
r2_res <- summary(modelo_baseball_res)$r.squared
n <- nobs(modelo_baseball_nores)
k <- 5 (Número de variables en el modelo no restringido)
q <- 3 (Número de variables que se quitaron)

# Calculo el estadistico F
f_stat_r2 <- ((r2_nores - r2_res) / (1 - r2_nores)) * ((n - k - 1) / q)
f_stat_r2

# Valores criticos
val_criticos <- c(qf(1 - 0.10, q, n - k - 1),
qf(1 - 0.05, q, n - k - 1),
qf(1 - 0.01, q, n - k - 1))
val_criticos

# Conclusion al 10%
if (f_stat > val_criticos[1]){
print("Se rechaza H0 en favor de H1 al 10% de significancia")
} else{
print("No se rechaza H0 al 10% de significancia")
}

# Conclusion al 5%
if (f_stat > val_criticos[2]){
print("Se rechaza H0 en favor de H1 al 5% de significancia")
} else{
print("No se rechaza H0 al 5% de significancia")
}

# Conclusion al 1%
if (f_stat > val_criticos[3]){
print("Se rechaza H0 en favor de H1 al 1% de significancia")
} else{
print("No se rechaza H0 al 1% de significancia")
}

# Calculo del valor p
p_val <- 1 - pf(f_stat, q, n - k - 1)
p_val

Escalado de datos

# Carga de base de datos
data(bwght)

# Definicion y seleccion de variables relevantes
base1 <- bwght %>%
mutate(peso_gr = bwght * 28.34952,
peso_kg = peso_gr / 1000,
cajetilla = cigs / 20) %>%
select(peso_gr, peso_kg, faminc, cigtax, cigs, cajetilla)

# Definicion modelos
modelo_peso_1 <- lm(formula = peso_gr ~ faminc + cigtax + cigs,
data = base1)
modelo_peso_2 <- lm(formula = peso_gr ~ faminc + cigtax + cajetilla,
data = base1)
modelo_peso_3 <- lm(formula = peso_kg ~ faminc + cigtax + cigs,
data = base1)
modelo_peso_4 <- lm(formula = peso_kg ~ faminc + cigtax + cajetilla,
data = base1)

# Presentacion de resultados
stargazer(modelo_peso_1, modelo_peso_2,
modelo_peso_3, modelo_peso_4,
digits = 3, type = "text", style = "qje",
report = c("vc*st"),
dep.var.labels = c("Peso en GR", "Peso en KG"),
covariate.labels = c("Ingreso Familiar",
"Impuestos",
"Cigarros",
"Cajetillas",
"Constante"))

Estandarización de datos

# Modelo de niveles
modelo_base <- lm(formula = price ~ crime + nox + rooms + dist + stratio,
data = base2)

# Modelo estandarizado
modelo_precio_casas <- lm(formula = scale(price) ~ 0 + scale(crime) + scale(nox) +
scale(rooms) + scale(dist) + scale(stratio),
data = base2)

# Presentacion resultados
stargazer(modelo_base, modelo_precio_casas, digits = 3,
type = "text", style = "io",
dep.var.labels = c("Precio (lvl)", "Precio (beta)"),
covariate.labels = c("Crimen", "Contaminación",
"Cuartos", "Distancia",
"Docente/Alumnado",
"Crimen STD", "Contaminación STD",
"Cuartos STD", "Distancia STD",
"Docente/Alumnado STD",
"Constante"))

Formas Logarítmicas

# Seleccion de variables relevantes
base3 <- hprice2 %>%
mutate(log_price = log(price),
log_nox = log(nox)) %>%
select(price, log_price, nox, log_nox)

# Modelos
modelo_lvl_lvl <- lm(formula = price ~ nox,
data = base3)
modelo_lvl_log <- lm(formula = price ~ log_nox,
data = base3)
modelo_log_lvl <- lm(formula = log_price ~ nox,
data = base3)
modelo_log_log <- lm(formula = log_price ~ log_nox,
data = base3)

# Presentacion resultados
stargazer(modelo_lvl_lvl, modelo_lvl_log,
modelo_log_lvl, modelo_log_log,
digits = 3, type = "text", style = "io",
dep.var.labels = c("Precio",
"LOG(Precio)"),
covariate.labels = c("Contaminación",
"LOG(Contaminación)",
"Constante"))

# Cuadricula 2x2
par(mfrow = c(2, 2))

# Grafica lvl-lvl
with(base3,
plot(nox, price, col = "dark gray",
main = "Level-level)",
xlab = "Contaminación",
ylab = "Precio"))
abline(modelo_lvl_lvl, col = "blue", lwd = 3)

# Grafica lvl-log
with(base3,
plot(log(nox), price, col = "dark gray",
main = "Level-log",
xlab = "LOG(Contaminación)",
ylab = "Precio"))
abline(modelo_lvl_log, col = "cyan", lwd = 3)

# Grafica log-lvl
with(base3,
plot(nox, log(price), col = "dark gray",
main = "Log-level",
xlab = "Contaminación",
ylab = "LOG(Precio)"))
abline(modelo_log_lvl, col = "red", lwd = 3)

# Grafica log-log
with(base3,
plot(log(nox), log(price), col = "dark gray",
main = "Log-log",
xlab = "LOG(Contaminación)",
ylab = "LOG(Precio)"))
abline(modelo_log_log, col = "orange", lwd = 3)

Formas Cuadráticas

# Selección de variables relevantes
base4 <- wage1 %>%
mutate(exper_sq = exper^2) %>%
select(wage, exper, exper_sq)

# Estimacion del modelo
modelo_experiencia <- lm(formula = wage ~ exper + exper_sq,
data = base4)

# Presentacion de resultados
stargazer(modelo_experiencia, digits = 3,
type = "text", style = "qje",
dep.var.labels = c("Salario"),
covariate.labels = c("Experiencia",
"Esperancia SQ",
"Constante"))
Los resultados obtenidos indican que existe un efecto marginal decreciente para la experiencia. Para poder ahondar más en estos retornos decrecientes necesitamos calcular el punto de quiebre (donde exper tiene efecto cero sobre el salario). Esto se hace de la siguiente manera:
x∗=β̂1/2β̂2=0.298/2(−0.006)=24.315
Esto indica que:

Quienes tienen menos de 24.315 años de experiencia reciben un efecto positivo por cada año adicional de experiencia.
Quienes tienen más de 24.315 años de experiencia reciben un efecto negativo por cada año adicional de experiencia.

R^2 para selección de modelos

# Selección de variables relevantes
base5 <- mlb1 %>%
mutate(log_salary = log(salary)) %>%
select(log_salary, years, games, runs, hits, hruns, rbis, bavg)

# Estimacion de los modelos
modelo_base1 <- lm(formula = log_salary ~ years + games + runs + hits + hruns + rbis + bavg,
data = base5)
modelo_base2 <- lm(formula = log_salary ~ games + runs + hits + bavg,
data = base5)
modelo_base3 <- lm(formula = log_salary ~ games + runs + bavg,
data = base5)

# Presentacion de resultados
stargazer(modelo_base1,
modelo_base2,
modelo_base3,
digits = 3,
type = "text", style = "qje",
dep.var.labels = c("Log(Salario)"),
covariate.labels = c("Antigüedad",
"Juegos",
"Carreras",
"Hits",
"Home Runs",
"Carreras Impulsadas",
"Promedio Bateo",
"Constante"))
De la tabla de resultados es importante notar que:

La R2
de los modelos es mayor en aquellos que incluyen más variables explicativas. Esto hace que el coeficiente de determinación sea una mala herramienta para decidir si una variable pertenece al modelo o no.
La R¯2
es una mejor herramienta (aunque no ideal) para elegir entre diversos modelos. En este caso el modelo 81) es el que tiene mayor capacidad para explicar el salario de los jugadores.
Es importante recordar que el criterio principal para decidir si una variable(s) pertenece(n) a un modelo es si significancia (conjunta) sobre la variable a explicar.

Variables Dummy para varias Categorías

# Seleccion de variables relevantes
base1 <- wage1 %>%
mutate(log_wage = log(wage)) %>%
select(log_wage, educ, female)

# Estimacion del modelo
discriminacion_simple <- lm(formula = log_wage ~ educ + female,
data = base1)

# Presentacion de resultados
stargazer(discriminacion_simple, digits = 3,
type = "text", style = "aer",
dep.var.labels = c("LOG(Salario)"),
covariate.labels = c("Escolaridad",
"Sexo", "Constante"))

# Grafica puntos
with(base1,
plot(educ, log_wage,
col = "dark gray",
pch = 16, cex = 0.50,
xlim = c(0, max(educ)),
ylim = c(0, 3.5),
main = "Discriminación Salarial por Género",
xlab = "Escolaridad",
ylab = "LOG(Salario)"))

# Recta Hombres
abline(a = coef(discriminacion_simple)[[1]],
b = coef(discriminacion_simple)[[2]],
col = "darkgreen", lwd = 3)

# Recta Mujeres
abline(a = coef(discriminacion_simple)[[1]] + coef(discriminacion_simple)[[3]],
b = coef(discriminacion_simple)[[2]],
col = "purple", lwd = 3)

# Legend
legend("topleft",
col = c("darkgreen", "purple"),
lwd = c(3, 3),
legend = c("Hombres", "Mujeres"))

base1m <- wage1 %>%
mutate(log_wage = log(wage),
exper_sq = exper^2,
tenure_sq = tenure^2) %>%
select(log_wage, educ, exper, exper_sq,
tenure, tenure_sq, female)

# Estimacion del modelo
discriminacion <- lm(formula = log_wage ~ educ + exper + exper_sq +
tenure + tenure_sq + female,
data = base1m)

# Presentacion de resultados
stargazer(discriminacion_simple,
discriminacion,
digits = 3,
type = "text", style = "aer",
dep.var.labels = c("LOG(Salario)"),
covariate.labels = c("Escolaridad",
"Experiencia",
"Experiencia SQ",
"Antigüedad",
"Antigüedad SQ",
"Sexo", "Constante"))

# Seleccion de variables relevantes
base2 <- beauty %>%
mutate(log_wage = log(wage),
exper_sq = exper^2,
looks_factor = as.factor(looks)) %>%
select(log_wage, educ, exper, exper_sq, looks_factor)

# Estimacion del modelo
atractive1 <- lm(formula = log_wage ~ educ + exper + exper_sq + looks_factor,
data = base2)

# Presentacion resultados
stargazer(atractive1, digits = 3,
type = "text", style = "aer",
dep.var.labels = c("LOG(Salario)"),
covariate.labels = c("Escolaridad", "Experiencia",
"Experiencia SQ", "Looks 2",
"Looks 3", "Looks 4", "Looks 5",
"Constante"))

# Redefinir grupo base
base2$looks_factor <- relevel(base2$looks_factor, "5")

# Estimacion del modelo
atractive5 <- lm(formula = log_wage ~ educ + exper + exper_sq + looks_factor,
data = base2)

# Presentacion resultados
stargazer(atractive5, digits = 3,
type = "text", style = "aer",
dep.var.labels = c("LOG(Salario)"),
covariate.labels = c("Escolaridad", "Experiencia",
"Experiencia SQ", "Looks 1",
"Looks 2", "Looks 3", "Looks 4",
"Constante"))

# Redefinir grupo base
base2$looks_factor <- relevel(base2$looks_factor, "3")

# Estimacion del modelo
atractive3 <- lm(formula = log_wage ~ educ + exper + exper_sq + looks_factor,
data = base2)

# Presentacion resultados
stargazer(atractive5, atractive1, atractive3,
digits = 3,
type = "text", style = "aer",
column.labels = c("Base Looks=5", "Base Looks=1", "Base Looks=3"),
dep.var.labels = c("LOG(Salario)"),
covariate.labels = c("Escolaridad", "Experiencia",
"Experiencia SQ", "Looks 1",
"Looks 2", "Looks 3",
"Looks 4", "Looks 5",
"Constante"))

# Seleccion de variables relevantes
base3 <- gpa3 %>%
filter(spring == 1) %>%
select(cumgpa, sat, hsperc, tothrs, female)

# Estimacion del modelo
modelo_diff <- lm(formula = cumgpa ~ female * (sat + hsperc + tothrs),
data = base3)

# Presentacion de resultados
stargazer(modelo_diff, digits = 5,
type = "text", style = "aer",
dep.var.labels = c("Promedio"),
covariate.labels = c("Mujer", "SAT", "Percentil Prepa",
"Horas Clase", "Mujer * SAT",
"Mujer * Percentil Prepa",
"Mujer * Horas Clase", "Constante"))

Heterocedasticidad

Prueba Breush-Pagan

# Estimacion del modelo
modelo_casas <- lm(formula = price ~ lotsize + sqrft + bdrms,
data = base1)

# Presentación de resultados
stargazer(modelo_casas, digits =3,
type = "text", style = "io",
dep.var.labels = c("Precio"),
covariate.labels = c("Tamaño Lote", "Tamaño Construcción",
"Número de Habitaciones", "Constante"))

# Paso 1:
# Recuperar los residuales al cuadrado
u2 <- resid(modelo_casas)^2

# Paso 2:
# Regresion auxiliar de los residuales
reg_auxiliar <- lm(formula = u2 ~ lotsize + sqrft + bdrms,
data = base1)

# Paso 3:
# Conclusion con estadistico F/valor-p
stargazer(reg_auxiliar, digits = 3,
type = "text", style = "io",
dep.var.labels = c("Residuales Cuadrados"),
covariate.labels = c("Tamaño Lote", "Tamaño Construcción",
"Número de Habitaciones", "Constante"))

Tenemos que el estadístico F
de la preba de significancia global de la regresión auxiliar es de 5.339 y su valor-p
es de 0.002.

Este valor del estadístico y del valor-p
de la prueba de significancia global nos llevan a rechazar la hipótesis nula, que en el caso de la prueba Breusch-Pagan es que el modelo cumple el supuesto RLM5.

Por tanto, al rechazar la hipótesis nula, estamos encontrando evidencia de que el modelo sufre de heterocedasticidad y, por tanto, se deberían tomar medidas correctivas, como el uso de errores estándar robustos.

Pruebas de White

# Paso 1:
# Recuperar los residuales al cuadrado
u2 <- resid(modelo_casas)^2

# Paso 2:
# Regresion auxiliar de los residuales
aux_white <- lm(formula = u2 ~ lotsize + sqrft + bdrms +
I(lotsize^2) + I(sqrft^2) + I(bdrms^2) +
lotsize * sqrft +
lotsize * bdrms +
sqrft * bdrms,
data = base1)

# Paso 3:
# Conclusion con estadistico F/valor-p
stargazer(aux_white, digits = 3,
type = "text", style = "io",
dep.var.labels = c("Residuales Cuadrados"),
covariate.labels = c("Tamaño Lote", "Tamaño Construcción",
"Número de Habitaciones",
"Tamaño Lote-SQ", "Tamaño Construcción-SQ",
"Número de Habitaciones-SQ",
"Lote * Construcción", "Lote * Habitaciones",
"Construcción * Habitaciones", "Constante"))

Tenemos que el estadístico F
de la preba de significancia global de la regresión auxiliar es de 5.387 y su valor-p
es de 0.00001013.

Este valor del estadístico y del valor-p
de la prueba de significancia global nos llevan a rechazar la hipótesis nula, que en el caso de la prueba de White es que el modelo cumple el supuesto RLM5.

Por tanto, al rechazar la hipótesis nula, estamos encontrando evidencia de que el modelo sufre de heterocedasticidad.

Notemos que esta prueba tiene algunos. Primero, requiere incorporar muchos regresores y es fácil equivocarse al omitir alguno de ellos. Segundo, al incorporar tantos regresores afecta los grados de libertad de la prueba F
.

Para resolver estos problemas, se cuenta con la prueba de White simplificada. Esta prueba también sirve para detectar heterocedasticidad pero usa los valores ajustados de la regresión en lugar de los regresores para encontrar relaciones entre los residuales y las variables explicativas.

Prueba White Simplificada

# Paso 1:
# Recuperar los residuales al cuadrado y valores ajustados
u2 <- resid(modelo_casas)^2
y_hat <- fitted(modelo_casas)

# Paso 2:
# Regresión auxiliar de residuales y valores ajustados
aux_whitesimp <- lm(formula = u2 ~ y_hat + I(y_hat^2))

# Paso 3:
# Conclusion con estadistico F/valor-p
stargazer(aux_whitesimp, digits = 3,
type = "text", style = "io",
dep.var.labels = c("Residuales Cuadrados"),
covariate.labels = c("Ajustados", "Ajustados-SQ",
"Constante"))

Tenemos que el estadístico F
de la preba de significancia global de la regresión auxiliar es de 9.639 y su valor-p
es de 0.0001687.

Este valor del estadístico y del valor-p
de la prueba de significancia global nos llevan a rechazar la hipótesis nula, que en el caso de la prueba de White simplificada es que el modelo cumple el supuesto RLM5.

Por tanto, al rechazar la hipótesis nula, estamos encontrando evidencia de que el modelo sufre de heterocedasticidad y, por tanto, debemos tomar medidas correctivas.

Errores estándar e inferencia Robusta a heterocedasticidad

# Presentación de resultados
stargazer(modelo_casas, digits =3,
type = "text", style = "io",
dep.var.labels = c("Precio"),
covariate.labels = c("Tamaño Lote", "Tamaño Construcción",
"Número de Habitaciones", "Constante"))

# Modelo base
(clasico <- coeftest(modelo_casas))

# Modelo BP
(bp <- coeftest(modelo_casas,
vcov = hccm(modelo_casas, type = "hc2")))

# Modelo White
(white <- coeftest(modelo_casas,
vcov = hccm(modelo_casas, type = "hc0")))

# Modelo White simplificado
(whitesimp <- coeftest(modelo_casas,
vcov = hccm))

# Presentacion de resultados
stargazer(clasico, bp, white, whitesimp,
digits = 5, type = "text", style = "io",
dep.var.labels = c("Precio"),
covariate.labels = c("Tamaño Lote",
"Tamaño Construcción",
"Número de Habitaciones",
"Constante"),
column.labels = c("Clásico", "Breusch-Pagan",
"White", "White Simplificado"))

Sin corrección por heterocedasticidad tenemos:

Coeficiente estimado: 0.00207
Error estándar: 0.00064
Estadístico t
: 3.2201
Valor-p
de la prueba: 0.00182
Conclusión: Se rechaza H0
Con corrección BP:

Coeficiente estimado: 0.00207
Error estándar: 0.00287
Estadístico t
: 0.71957
Valor-p
de la prueba: 0.47378
Conclusión: No se rechaza H0
Con corrección de White:

Coeficiente estimado: 0.00207
Error estándar: 0.00122
Estadístico t
: 1.69117
Valor-p
de la prueba: 0.09451
Conclusión: No se rechaza H0
Con corrección de White simplificada:

Coeficiente estimado: 0.00207
Error estándar: 0.00715
Estadístico t
: 0.28925
Valor-p
de la prueba: 0.7731
Conclusión: No se rechaza H0

# Carga de datos y eliminación de datos faltantes
data(mroz)
mroz %<>%
filter(!is.na(wage)) %>%
select(wage, educ, fatheduc, motheduc)

# Cálculo manual de estimadores de MCO y de IV
educ_mco <- with(mroz, cov(log(wage), educ) / var(educ))
educ_iv <- with(mroz, cov(log(wage), fatheduc) / cov(educ, fatheduc))
c(MCO = educ_mco,
IV =educ_iv)

# Regresión de MCO usando lm
mco_reg <- lm(formula = log(wage) ~ educ,
data = mroz)
# Regresión de IV usando ivreg
iv_reg <- ivreg(formula = log(wage) ~ educ | fatheduc,
data = mroz)

# Presentación de resultados usando stargazer
stargazer(mco_reg, iv_reg, type = "text",
digits = 3, style = "io",
model.names = FALSE,
dep.var.labels = c("Log(wage)"),
column.labels = c("OLS", "IV"),
covariate.lables = c("Educación","Constante"))

# Regresión de IV con motheduc como instrumento
iv_reg2 <- ivreg(formula = log(wage) ~ educ | motheduc,
data = mroz)
# Presentación de los dos modelos de IV
stargazer(iv_reg, iv_reg2, type = "text",
digits = 3, style = "io",
model.names = FALSE,
dep.var.labels = c("Log(wage)"),
column.labels = c("Padre", "Madre"),
covariate.labels = c("Educación", "Constante"))

# Carga de datos
data(card)
card %<>%
select(wage, educ, nearc4, exper, black, smsa,
south, smsa66, reg662, reg663, reg664,
reg665, reg666, reg667, reg668, reg669)

# Estimación de la forma reducida del modelo
mod_red <- lm(educ ~ nearc4 + exper + I(exper^2) + black +
smsa + south + smsa66 + reg662 + reg663 +
reg664 + reg665 + reg666 + reg667 +
reg668 + reg669,
data = card)
# Prueba t sobre el coeficiente del instrumento
linearHypothesis(mod_red, c("nearc4 = 0"))

# Estimación del modelo por MCO
mod_lm <- lm(log(wage) ~ educ + exper + I(exper^2) +
black + smsa + south + smsa66 +
reg662 + reg663 + reg664 + reg665 +
reg666 + reg667 + reg668 + reg669,
data = card)
# Estimación del modelo pot IV
mod_iv <- ivreg(log(wage) ~ educ + exper + I(exper^2) +
black + smsa + south + smsa66 +
reg662 + reg663 + reg664 + reg665 +
reg666 + reg667 + reg668 + reg669 |
nearc4 + exper + I(exper^2) +
black + smsa + south + smsa66 +
reg662 + reg663 + reg664 + reg665 +
reg666 + reg667 + reg668 + reg669,
data = card)

# Presentación de los resultados
stargazer(mod_lm, mod_iv, type = "text",
digits = 3, style = "io",
model.names = FALSE,
dep.var.labels = c("log(wage)"),
column.labels = c("OLS", "IV"),
keep = c("educ", "exper", "black"),
covariate.labels = c("Educación", "Experiencia",
"Afrodescendiente", "Constante"))

Minimos Cuadrados de 2 Etapas

# Carga de datos
data(mroz)
mroz %<>%
filter(!is.na(wage)) %>%
select(wage, educ, fatheduc, motheduc, exper)
# Modelo OLS
mod_lm <- lm(log(wage) ~ educ + exper + I(exper^2),
data = mroz)
# Modelo IV con fatheduc como instrumento
mod_iv1 <- ivreg(log(wage) ~ educ + exper + I(exper^2) |
fatheduc + exper + I(exper^2),
data = mroz)
# Modelo IV con motheduc como instrumento
mod_iv2 <- ivreg(log(wage) ~ educ + exper + I(exper^2) |
motheduc + exper + I(exper^2),
data = mroz)
# Presentación de los resultados
stargazer(mod_lm, mod_iv1, mod_iv2, type = "text",
digits = 3, style = "io",
model.names = FALSE,
column.labels = c("OLS", "Father IV", "Mother IV"),
dep.var.labels = c("log(wage)"),
covariate.labels = c("Educación", "Experiencia",
"Experiencia-SQR", "Constante"))

# Paso 1: Forma reducida y prueba de relevancia
mod_red <- lm(educ ~ exper + I(exper^2) + fatheduc + motheduc,
data = mroz)
linearHypothesis(mod_red, c("fatheduc = 0", "motheduc = 0"))

# Paso 2: Forma estructural
mod_est <- lm(log(wage) ~ fitted(mod_red) + exper + I(exper^2),
data = mroz)

mod_iv <- ivreg(log(wage) ~ educ + exper + I(exper^2) |
fatheduc + motheduc + exper + I(exper^2),
data = mroz)

# Presentacion de los resultados
stargazer(mod_est, mod_iv, type = "text",
digits = 3, style = "io",
model.names = FALSE,
column.labels = c("MANUAL", "IVREG"),
dep.var.labels = c("Log(wage)"),
covariate.labels = c("Educación (Fitted)",
"Educación", "Experiencia",
"Experiencia-SQR", "Constante"))

Pruebas de Endogeneidad y sobreidentificación

# Carga de los datos
data(mroz)
mroz %<>%
filter(!is.na(wage)) %>%
select(wage, educ, fatheduc, motheduc, exper)
# Paso 1: Forma reducida
mod_red <- lm(educ ~ exper + I(exper^2) + fatheduc + motheduc,
data = mroz)
# Paso 2: Forma estructural con residuales
mod_est <- lm(log(wage) ~ educ + exper + I(exper^2) +
resid(mod_red),
data = mroz)
# Se presentan los resultados
stargazer(mod_est, type = "text",
digits = 3, style = "io",
dep.var.labels = c("Log(wage)"),
covariate.labels = c("Educación", "Experiencia",
"Experiencia-SQR",
"Residuales (forma reducida)",
"Constante"))

# Se realiza la prueba t
linearHypothesis(mod_est, c("resid(mod_red) = 0"))

# Paso 1: Estimación de 2SLS y residuales
mod_2sls <- ivreg(log(wage) ~ educ + exper +I(exper^2) |
fatheduc + motheduc + exper + I(exper^2),
data = mroz)
r_2sls <- resid(mod_2sls)
# Paso 2: Regresión sobre los residuales
mod_resid <- lm(r_2sls ~ exper + I(exper^2) +
fatheduc + motheduc,
data = mroz)
# Paso 3: Cálculo del estadístico y conclusión
r_sq <- summary(mod_resid)$r.squared
n <- nobs(mod_resid)
(test_stat <- n * r_sq)

q <- 2 - 1
(p_val <- 1 - pchisq(test_stat, q))