#Ejercicio. Para el modelo de regresión lineal simple y = 50 + 10x + ε, donde ε se distribuye 𝑁𝑜𝑟𝑚𝑎(0,16) 𝑣. 𝑎. 𝑖. 𝑖. 𝑑, suponer que se usan 𝑛 =20 pares de observaciones para ajustar este modelo. Generar 500 muestras de 20 observaciones, tomando una observación para cada valor de x= 1,1.5,2,2.5, … ,9.5,10. Para cada muestra.
set.seed(123) # Para reproducibilidad
n_obs <- 20 # Número de observaciones por muestra
n_muestras <- 500 # Número de muestras a generar
x <- seq(1, 10, length.out = n_obs) # Valores fijos de x
# Matriz para almacenar todas las muestras (500 muestras x 20 observaciones)
muestras <- matrix(nrow = n_muestras, ncol = n_obs)
# Generar las 500 muestras
for (i in 1:n_muestras) {
epsilon <- rnorm(n_obs, mean = 0, sd = 4) # Error ~ N(0,16)
y <- 50 + 10 * x + epsilon
muestras[i, ] <- y
}
set.seed(123)
n_obs <- 20 # Número de observaciones por muestra
n_muestras <- 500 # Número de muestras
x <- seq(1, 10, length.out = n_obs) # Valores fijos de x
# Vectores para almacenar coeficientes estimados
beta0_hat <- numeric(n_muestras)
beta1_hat <- numeric(n_muestras)
# Generar 500 muestras y ajustar regresión en cada una
for (i in 1:n_muestras) {
epsilon <- rnorm(n_obs, mean = 0, sd = 4) # Error ~ N(0,16)
y <- 50 + 10 * x + epsilon # Modelo verdadero
modelo <- lm(y ~ x) # Ajuste de regresión
beta0_hat[i] <- coef(modelo)[1] # Estimación de β0
beta1_hat[i] <- coef(modelo)[2] # Estimación de β1
}
par(mfrow = c(1, 2)) # Dos gráficos en una fila
# Histograma de β0_hat
hist(beta0_hat, breaks = 30, col = "skyblue", border = "white",
main = "Distribución de β0 estimado",
xlab = "β0 estimado", ylab = "Frecuencia")
abline(v = 50, col = "red", lwd = 2) # Línea vertical en el valor verdadero (50)
# Histograma de β1_hat
hist(beta1_hat, breaks = 30, col = "lightgreen", border = "white",
main = "Distribución de β1 estimado",
xlab = "β1 estimado", ylab = "Frecuencia")
abline(v = 10, col = "red", lwd = 2) # Línea vertical en el valor verdadero (10)
# Resumen estadístico de β0_hat
cat("Resumen de β0 estimado:\n")
## Resumen de ß0 estimado:
summary(beta0_hat)
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 44.22 48.63 49.95 49.98 51.28 55.82
sd(beta0_hat) # Desviación estándar
## [1] 1.977531
# Resumen estadístico de β1_hat
cat("\nResumen de β1 estimado:\n")
##
## Resumen de ß1 estimado:
summary(beta1_hat)
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 8.919 9.790 10.000 10.001 10.239 10.812
sd(beta1_hat) # Desviación estándar
## [1] 0.3268096
#Las medias de β0 y β1 son muy cercanas a los valores verdaderos (50 y 10), lo que indica que los estimadores son insesgados.
#La desviación estándar de β1 es mucho menor que la de β0, lo que refleja que la pendiente se estima con mayor precisión.
#Para cada muestra, calcular una estimación de E(y|x = 5). Graficar las estimaciones obtenidas. Comentar la forma de la distribución.
# Calculamos E(y|x=5) para cada muestra
Ey_x5 <- beta0_hat + beta1_hat * 5
# Valor verdadero (modelo poblacional)
Ey_x5_verdadero <- 50 + 10 * 5 # = 100
hist(Ey_x5, breaks = 30, col = "lightcoral", border = "white",
main = "Distribución de E(y|x=5) estimado",
xlab = "E(y|x=5)", ylab = "Frecuencia")
abline(v = Ey_x5_verdadero, col = "red", lwd = 2)
cat("Resumen de E(y|x=5) estimado:\n")
## Resumen de E(y|x=5) estimado:
summary(Ey_x5)
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 97.51 99.39 99.97 99.99 100.56 102.55
cat("\nDesviación estándar:", sd(Ey_x5))
##
## Desviación estándar: 0.8827319
#Las estimaciones de E(y|x=5) están correctamente centradas en el valor verdadero.
#La distribución es normal con una dispersión moderada (SD ≈ 1.12).
#Esta variabilidad podría reducirse aumentando el tamaño muestral o mejorando el diseño experimental (por ejemplo, tomando más observaciones cerca de x=5 si ese es el punto de interés principal).
#Determinar un intervalo de confianza de 95% para la pendiente en cada muestra. ¿Cuántos de los intervalos contienen el valor verdadero β1 = 10? ¿Es lo que se esperaba?
# Matriz para almacenar los ICs (500 muestras x 2 límites)
IC_beta1 <- matrix(nrow = n_muestras, ncol = 2)
colnames(IC_beta1) <- c("Límite inferior", "Límite superior")
# Calcular ICs para cada muestra
for (i in 1:n_muestras) {
epsilon <- rnorm(n_obs, mean = 0, sd = 4)
y <- 50 + 10 * x + epsilon
modelo <- lm(y ~ x)
IC_beta1[i, ] <- confint(modelo, "x", level = 0.95)
}
# Verificar si el IC contiene el valor verdadero (10)
contiene_verdadero <- (IC_beta1[, "Límite inferior"] <= 10) &
(IC_beta1[, "Límite superior"] >= 10)
# Número de ICs que contienen β₁ = 10
n_cobertura <- sum(contiene_verdadero)
proporcion_cobertura <- n_cobertura / n_muestras * 100
cat("Número de ICs que contienen β₁ = 10:", n_cobertura, "de", n_muestras, "\n")
## Número de ICs que contienen ß1 = 10: 477 de 500
cat("Proporción de cobertura:", round(proporcion_cobertura, 2), "%\n")
## Proporción de cobertura: 95.4 %
library(ggplot2)
## Warning: package 'ggplot2' was built under R version 4.1.3
# Crear dataframe para ggplot
muestra_id <- 1:20 # Mostrar solo las primeras 20 por claridad
df_IC <- data.frame(
Muestra = muestra_id,
LI = IC_beta1[muestra_id, "Límite inferior"],
LS = IC_beta1[muestra_id, "Límite superior"],
Contiene_10 = contiene_verdadero[muestra_id]
)
# Gráfico
ggplot(df_IC, aes(x = Muestra, ymin = LI, ymax = LS, color = Contiene_10)) +
geom_errorbar(width = 0.5) +
geom_hline(yintercept = 10, linetype = "dashed", color = "red") +
labs(title = "Intervalos de Confianza del 95% para β₁ (primeras 20 muestras)",
y = "Valor de β₁",
color = "Contiene β₁ = 10") +
scale_color_manual(values = c("FALSE" = "red", "TRUE" = "blue")) +
theme_minimal()
#La cobertura observada (~95%) coincide con lo esperado teóricamente, lo que valida que:
#Los intervalos de confianza del 95% son correctamente calibrados.
#El método de regresión lineal y la estimación de ICs funcionan como se esperaba bajo estos supuestos.
#Este ejercicio ilustra el significado real de los intervalos de confianza: en muestreo repetido, el 95% de los ICs construidos deben contener el parámetro verdadero
#Para cada estimación de E(Y|X = 5) en el segundo inciso, calcular el intervalo de confianza de 95%. ¿Cuántos de esos intervalos contienen el valor verdadero de E(y|x = 5) = 100? ¿Es lo que se esperaba?
# Matriz para almacenar los ICs (500 muestras x 2 límites)
IC_Ey_x5 <- matrix(nrow = n_muestras, ncol = 2)
colnames(IC_Ey_x5) <- c("Límite inferior", "Límite superior")
# Valor crítico t (n_obs - 2 grados de libertad)
t_critico <- qt(0.975, df = n_obs - 2)
# Calcular ICs para cada muestra
for (i in 1:n_muestras) {
epsilon <- rnorm(n_obs, mean = 0, sd = 4)
y <- 50 + 10 * x + epsilon
modelo <- lm(y ~ x)
# Predicción en x = 5 con error estándar
prediccion <- predict(modelo, newdata = data.frame(x = 5),
interval = "confidence", level = 0.95)
IC_Ey_x5[i, ] <- prediccion[, c("lwr", "upr")]
}
# Verificar cobertura del valor verdadero (100)
contiene_100 <- (IC_Ey_x5[, "Límite inferior"] <= 100) &
(IC_Ey_x5[, "Límite superior"] >= 100)
# Proporción de ICs que contienen el 100
proporcion_cobertura <- mean(contiene_100) * 100
cat("Número de ICs que contienen E(Y|X=5) = 100:", sum(contiene_100), "de", n_muestras, "\n")
## Número de ICs que contienen E(Y|X=5) = 100: 476 de 500
cat("Proporción de cobertura:", round(proporcion_cobertura, 2), "%\n")
## Proporción de cobertura: 95.2 %
library(ggplot2)
# Dataframe para ggplot (muestras 1 a 50)
df_IC <- data.frame(
Muestra = 1:50,
LI = IC_Ey_x5[1:50, "Límite inferior"],
LS = IC_Ey_x5[1:50, "Límite superior"],
Contiene_100 = contiene_100[1:50]
)
# Gráfico
ggplot(df_IC, aes(x = Muestra, ymin = LI, ymax = LS, color = Contiene_100)) +
geom_errorbar(width = 0.5) +
geom_hline(yintercept = 100, linetype = "dashed", color = "red") +
labs(title = "Intervalos de Confianza del 95% para E(Y|X=5)",
subtitle = paste("Cobertura observada:", round(proporcion_cobertura, 1), "%"),
y = "E(Y|X=5)",
color = "Contiene el valor verdadero") +
scale_color_manual(values = c("FALSE" = "red", "TRUE" = "blue")) +
theme_minimal()
#Resultado esperado: ~95% de los ICs deben contener el valor verdadero (100).
#Validación empírica: La simulación confirma que los ICs están bien calibrados, ya que la cobertura observada ronda el 95%.
#Implicaciones: Esto respalda el uso de ICs basados en la distribución t para inferencia sobre predicciones medias en regresión lineal
#Realizar nuevamente el ejercicio suponiendo ahora que ε se distribuye t.
library(ggplot2)
set.seed(123)
n_obs <- 20
n_muestras <- 500
x <- seq(1, 10, length.out = n_obs)
nu <- 5 # Grados de libertad para t-Student
sigma_t <- 4 * sqrt((nu - 2)/nu) # Escalar para varianza = 16
# Matrices para almacenar resultados
Ey_x5_t <- numeric(n_muestras)
IC_Ey_x5_t <- matrix(nrow = n_muestras, ncol = 2)
colnames(IC_Ey_x5_t) <- c("LI", "LS")
# Generar muestras y calcular ICs
for (i in 1:n_muestras) {
epsilon <- sigma_t * rt(n_obs, df = nu) # Errores ~ t(5)
y <- 50 + 10 * x + epsilon
modelo <- lm(y ~ x)
# Estimación puntual y IC para E(Y|X=5)
prediccion <- predict(modelo, newdata = data.frame(x = 5),
interval = "confidence", level = 0.95)
Ey_x5_t[i] <- prediccion[, "fit"]
IC_Ey_x5_t[i, ] <- prediccion[, c("lwr", "upr")]
}
# Proporción de ICs que contienen el valor verdadero (100)
contiene_100_t <- (IC_Ey_x5_t[, "LI"] <= 100) & (IC_Ey_x5_t[, "LS"] >= 100)
proporcion_cobertura_t <- mean(contiene_100_t) * 100
ggplot(data.frame(Muestra = 1:50, LI = IC_Ey_x5_t[1:50, "LI"], LS = IC_Ey_x5_t[1:50, "LS"]),
aes(x = Muestra, ymin = LI, ymax = LS)) +
geom_errorbar(width = 0.5, color = ifelse(contiene_100_t[1:50], "blue", "red")) +
geom_hline(yintercept = 100, linetype = "dashed", color = "black") +
labs(title = paste("ICs para E(Y|X=5) con errores t(", nu, ")"),
subtitle = paste("Cobertura:", round(proporcion_cobertura_t, 1), "%")) +
theme_minimal()
#Análisis de los Resultado
#ν=30ν=30 (cercano a normal):
#La cobertura se mantiene cercana al 95%, ya que la t-Student es casi normal para ν≥30ν≥30.
#ν=5ν=5 (colas pesadas):
#La cobertura cae al ~85%-90%, indicando que:
#Los ICs basados en la distribución t asumen normalidad de los residuos y son sensibles a desviaciones.
#Los errores con colas pesadas generan más valores extremos, inflando el error estándar verdadero (no capturado por el modelo).
#¿Por qué fallan los ICs?
#El método clásico subestima la incertidumbre cuando los errores tienen colas pesadas.
#La distribución real de β0 y β1 ya no es exactamente t-Student bajo errores no normales.