# ======================================================== #
#     Tema 9: Problemas comunes en modelos de regresión    #
# ======================================================== #

# ========================
# CARGA DE DATOS ---------
# ========================
library(readxl)
base_enemdu <- read_excel("G:/Mi unidad/3. Trabajos/8. Senamhi  Análisis de datos en R/4. Curso R/Módulo 2/Ejercicios 6/base_simulada_enemdu.xlsx")


# ==========================
# MODELO BASE -------------
# ==========================

base_enemdu <- base_enemdu %>%
  mutate(ln_ingreso = log(ingreso))

modelo_1 <- lm(ln_ingreso ~ sexo + horas + experiencia + edad, data = base_enemdu)
summary(modelo_1)

modelo_2 <- lm(ln_ingreso ~ sexo + log(horas) + log(experiencia) + log(edad), data = base_enemdu)
summary(modelo_2)

# **************************
# 1. MULTICOLINEALIDAD -----
# **************************

# Matriz de correlación
cor(base_enemdu)

# VIF
install.packages("car")
library(car)

vif(model_1)

# Caso simulado con multicolinealidad
set.seed(123)
n <- 100

x1 <- rnorm(n)
x2 <- 0.9*x1 + rnorm(n, sd=0.2)  # alta correlación
x3 <- rnorm(n)

y <- 2 + 3*x1 + 2*x2 + 1*x3 + rnorm(n)

data_multi <- data.frame(y, x1, x2, x3)

cor(data_multi)

model_2 <- lm(y ~ x1 + x2 + x3, data=data_multi)
summary(model_2)
vif(model_2)


# ***************************
# 2. HETEROCEDASTICIDAD -----
# ***************************

install.packages("lmtest")
library(lmtest)

# Test Breusch-Pagan
bptest(model_1)

# Gráfico de residuos
plot(data_wage$exper, resid(model_1))

install.packages("sandwich")
library(sandwich)

# Errores robustos
coeftest(model_1, vcov = vcovHC(model_1))

# Caso simulado
set.seed(123)
n <- 100

x <- runif(n, 1, 10)
error <- rnorm(n, sd = x)
y <- 2 + 3*x + error

data_het <- data.frame(y, x)

model_het <- lm(y ~ x, data=data_het)
summary(model_het)

bptest(model_het)
bptest(model_het, ~ x + I(x^2))  # White aproximado

coeftest(model_het, vcov = vcovHC(model_het, type="HC1"))


# ***************************
# 3. AUTOCORRELACIÓN -------
# ***************************

# Test Durbin-Watson
dwtest(model_1)

# Simulación de autocorrelación
set.seed(123)
n <- 100

u <- arima.sim(model = list(ar=0.8), n=n)
x <- rnorm(n)
y <- 1 + 2*x + u

data_auto <- data.frame(y, x)

model_auto <- lm(y ~ x, data=data_auto)

# Test de Durbin–Watson
dwtest(model_auto)
bgtest(model_auto, order = 2)

# Extraer residuos
residuos <- resid(model_auto)

# Gráfico en el tiempo
plot(residuos, type = "l",
     main = "Residuos en el tiempo",
     xlab = "Tiempo",
     ylab = "Residuos")

# Línea horizontal en 0
abline(h = 0, col = "red")

# ***************************
# 4. ENDOGENEIDAD ----------
# ***************************

# Simulación
set.seed(123)
n <- 100

z <- rnorm(n)  # instrumento
u <- rnorm(n)

x <- 0.5*z + u
y <- 2 + 3*x + u

data_endo <- data.frame(y, x, z)

model_endo <- lm(y ~ x, data=data_endo)
summary(model_endo)

# Variables instrumentales
install.packages("AER")
library(AER)

iv_model <- ivreg(y ~ x | z, data=data_endo)
summary(iv_model)


# ***************************
# 5. VARIABLES OMITIDAS -----
# ***************************

# Modelo mal especificado
model_omit <- lm(lwage ~ educ + exper, data=data_wage)
summary(model_omit)

# Comparación con modelo completo
summary(model_1)


# ***************************
# 6. EFECTOS NO OBSERVADOS -
# ***************************

# Simulación panel simple
set.seed(123)
n <- 50
t <- 5

id <- rep(1:n, each=t)
tiempo <- rep(1:t, n)

alpha <- rnorm(n)  # efecto no observado
alpha_rep <- rep(alpha, each=t)

x <- rnorm(n*t)
u <- rnorm(n*t)

y <- 1 + 2*x + alpha_rep + u

data_panel <- data.frame(id, tiempo, y, x)

install.packages("plm")
library(plm)

# Modelo pooled
pool <- plm(y ~ x, data=data_panel, model="pooling")

# Modelo efectos fijos
fe <- plm(y ~ x, data=data_panel, model="within")

summary(pool)
summary(fe)


# ***************************
# 7. OUTLIERS --------------
# ***************************

# Distancia de Cook
cooksd <- cooks.distance(model_1)

# Visualización
plot(cooksd, type="h")

# Observaciones influyentes
which(cooksd > (4/length(cooksd)))


# ***************************
# 8. NO SIGNIFICANCIA -------
# ***************************

summary(modelo_1)

# Revisar p-values directamente
coef(summary(modelo_1))

# Intervalos de confianza
confint(modelo_1)


# ==========================
# FIN DEL SCRIPT ----------
# ==========================