Planteamiento

Suponer que se conocen los parámetros verdaderos de la distribución de \(Y\) condicional a \(X\), de manera que \[ E(Y \,\vert\, X = x) = 30 - 2x \qquad \text{y} \qquad V(Y \,\vert\, X = x) = 4 \]

Se genera una muestra tamaño \(n = 30\) de \(X\) que se mantendrá fija en todo el ejemplo

x <- rnorm(30, 10, 3); x
##  [1]  7.143442  9.725740  9.399617 11.402659  6.294711  9.770939 13.138348
##  [8]  7.332940 12.209379  7.579736  9.122284 15.045028  7.115359  8.232422
## [15]  9.703490  9.397131 11.568066 11.540788 13.463723 13.251499 10.068707
## [22] 12.741647  9.018621  9.939268  8.893413  7.341928 12.928035  4.413202
## [29]  7.695336 11.915824

Para cada valor de \(X\) se genera una \(Y\) y se ajusta el modelo que servira como base para las inferencias posteriores

y <- 30 - 2*x + rnorm(30, 0, 2)
modelo <- lm(y ~ x)
summary(modelo)
## 
## Call:
## lm(formula = y ~ x)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -3.7317 -1.4366 -0.1172  1.8448  4.3523 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  30.7840     1.7403   17.69  < 2e-16 ***
## x            -2.0471     0.1704  -12.02 1.44e-12 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 2.304 on 28 degrees of freedom
## Multiple R-squared:  0.8376, Adjusted R-squared:  0.8318 
## F-statistic: 144.4 on 1 and 28 DF,  p-value: 1.441e-12

En la siguiente gráfica se representan los datos y las rectas de regresión ajustada y verdadera.

Intervalo de confianza para \(E(Y\,\vert\, X=x)\)

Suponer que se tiene interés en hacer inferencias sobre el valor esperado de \(Y\) dado \(X = 14\) (este valor es arbitrario, pero está en el rango válido para el modelo). Se puede utilizar la función predict para realizar tales inferencias. El intervalo que se reporta es de confianza 95%.

predict(modelo, data.frame(x = 14), interval = 'confidence')
##        fit       lwr      upr
## 1 2.124576 0.4583902 3.790762

Como se conoce el modelo verdadero, se sabe que \(\mu_{14} = E(Y\,\vert\, 14) = 2\).

mu_14 <- 30 - 2 * 14; mu_14
## [1] 2

Se puede observar que con esta muestra, el intervalo de confianza calculado sí contiene al valor verdadero de \(\mu_{14}\).

Con el siguiente código se evalua la cobertura del intervalo de confianza 95% para \(\mu_{14}\), a partir de simular \(m = 5,000\) muestras de \(Y\). Se debe observar que no se generan nuevos valores de \(X\) y que solamente se generan valores de \(Y\).

cobertura <- c()
for(k in 1:5000)
{
  y_sim <- 30 - 2*x + rnorm(30, 0, 2)
  modelo_sim <- lm(y_sim ~ x)
  ic_sim <- predict(modelo_sim, data.frame(x = 14), interval = 'confidence')[2:3]
  cobertura[k] <- (ic_sim[1] <= mu_14) & (mu_14 <= ic_sim[2])
}
confianza <- 100 * sum(cobertura) / 5000
confianza
## [1] 95.44

Se puede ver cómo el porcentaje de cobertura con las \(5,000\) muestras es cercano al 95% especificado en la construcción del intervalo de confianza.

Intervalo de predicción para \(Y \,\vert\, X = x\).

Suponer que ahora se tiene interés se tiene interés es hacer una predicción para una nueva observación de \(Y\) dado \(X=14\). De nuevo se puede utilizar la función predict para simplificar las cosas. El intervalo que se calcula es predicción de confianza 95%.

predict(modelo, data.frame(x = 14), interval = 'prediction')
##        fit       lwr      upr
## 1 2.124576 -2.880019 7.129172

Como se trata de un ejemplo de simulación, podemos generar una nueva observación y comprobar si está incluída en el intervalo anterior.

nueva_obs <- function(x) 30 - 2*x + rnorm(1, 0, 2)
nueva_obs(14)
## [1] -0.06025884

Se puede observar que con la muestra inicial el intervalo de predicción de confianza 95% para una nueva observación de \(Y\) dado \(X = 14\), sí contiene a la observación generada.

Con el siguiente código se evalua la cobertura del intervalo de confianza 95% para una nueva observación de \(Y\) dado \(X = 14\), a partir de simular \(m = 5,000\) muestras de \(Y\) y en cada caso una nueva observación. De nuevo, es importante notar que no se generan nuevos valores de \(X\) y que solamente se generan valores de \(Y\).

cobertura <- c()
for(k in 1:5000)
{
  y_sim <- 30 - 2*x + rnorm(30, 0, 2)
  nueva_sim <- nueva_obs(14)
  modelo_sim <- lm(y_sim ~ x)
  ip_sim <- predict(modelo_sim, data.frame(x = 14), interval = 'prediction')[2:3]
  cobertura[k] <- (ip_sim[1] <= nueva_sim) & (nueva_sim <= ip_sim[2])
}
confianza <- 100 * sum(cobertura) / 5000
confianza
## [1] 94.88

Se puede ver cómo el porcentaje de cobertura con las \(5,000\) muestras es cercano al 95% especificado en la construcción del intervalo de predicción.

Intervalos simultáneos

Suponer que se tiene interés en construir intervalos simultáneos de confianza 95% para el valor esperado de \(Y\) dado \(X = x\), para \(x \in \{5, 7, 9, 11, 13\}\). Si se utiliza el método de Bonferroni, las inferencias se puede hacer fácilmente usando la función predict, solamente hace falta hacer un ajuste en el nivel de significancia.

Recordando, el método de Bonferroni para construir \(k\) intervalos simultáneos de confianza \(100(1-\alpha)\%\) indica que se debe usar el cuantil \(1-\alpha/2k\) de la distribución \(t_{(n-2)}\), pero esto es equivalente a constuir los \(k\) intervalos individuales de confianza \(1-\alpha/k\). Para este ejemplo, para construir 5 intervalos simultáneos de confianza 95%, se debe construir cada uno de confianza 99%.

predict(modelo, data.frame(x = seq(5, 13, 2)), interval = 'confidence', level = 0.99)
##         fit       lwr       upr
## 1 20.548491 17.960120 23.136861
## 2 16.454288 14.656711 18.251864
## 3 12.360084 11.120891 13.599278
## 4  8.265881  6.995993  9.535770
## 5  4.171678  2.310945  6.032411

Como se conocen los verdaderos parámetros del modelo es posible saber si cada intervalo contiene al parámetro para el que fue calculado.

m_x <- 30 - 2 * seq(5, 13, 2); m_x
## [1] 20 16 12  8  4

Se puede observar como los 5 intervalos contienen cada uno al parámetro para el que fueron calculados, esto significa que la muestra original fue una buena muestra.

Ahora se evaluará con simulación la confianza real de los intervalos simultáneos obtenidos por el método de Bonferroni, además se comparará la confianza conjunta en el caso en que los intervalos sean individuales. Nuevamente se consideran \(m = 5,000\) muestras para la simulación.

cobertura_ind <- c()
cobertura_sim <- c()
for(k in 1:5000)
{
  y_sim <- 30 - 2*x + rnorm(30, 0, 2)
  modelo_sim <- lm(y_sim ~ x)
  ic_ind <- predict(modelo_sim, data.frame(x = seq(5, 13, 2)), interval = 'confidence', level = 0.95)
  ic_sim <- predict(modelo_sim, data.frame(x = seq(5, 13, 2)), interval = 'confidence', level = 0.99)
  
  cobertura_ind[k] <- all(apply(cbind(ic_ind, m_x), 1, function(x) (x[2] <= x[4]) & (x[4] <= x[3])))
  cobertura_sim[k] <- all(apply(cbind(ic_sim, m_x), 1, function(x) (x[2] <= x[4]) & (x[4] <= x[3])))
}
confianza_ind <- 100 * sum(cobertura_ind) / 5000; confianza_ind
## [1] 86.62
confianza_sim <- 100 * sum(cobertura_sim) / 5000; confianza_sim
## [1] 96.94

Los resultados muestran que la confianza conjunta de los intervalos individuales con una confianza incorrecta, es mucho menor al 95% con el que fueron construidos y por otro lado, la confianza conjunta de los intervalos con la confianza corregida según el método de Bonferroni es mayor al 95% requerido, lo cual está bien, puesto que para específicar que la confianza es 95% se debe cumplir que la cobertura sea al menos ese porcentaje.

Bandas de confianza para la recta de regresión

Por último se evaluará la cobertura de las bandas de confianza, aunque con menos detalle que en los casos anteriores. Primero se calculan las cantidades fijas (que sólo dependen de \(X\)).

n <- length(x)
xbar <- mean(x)
sxx <- sum((x - xbar)**2)

Posteriormente se definen las funciones que calculan los límites de la bandas de confianza. Las funciones tienen tres argumentos: el valor de \(x\), los parámetros del modelo ajustado (esto cambia con cada muestra) y el nivel de confianza.

l_x <- function(x_n, ajuste, alfa)
{
  b0 <- coef(ajuste)[1]
  b1 <- coef(ajuste)[2]
  s2 <- summary(ajuste)$sigma**2
  (b0 + b1*x_n) - 2*qf(1 - alfa, 2, n-2)*sqrt(1/n + (x_n - xbar)^2/sxx)*sqrt(s2)
}

u_x <- function(x_n, ajuste, alfa)
{
  b0 <- coef(ajuste)[1]
  b1 <- coef(ajuste)[2]
  s2 <- summary(ajuste)$sigma**2
  (b0 + b1*x_n) + 2*qf(1 - alfa, 2, n-2)*sqrt(1/n + (x_n - xbar)^2/sxx)*sqrt(s2)
}

En la siguiente gráfica se muestran la recta de regresión verdadera, la recta de regresión ajustada con la muestra inicial y 10 bandas de confianza 75% para la recta de regresión veradera. Se elige \(\alpha = 0.25\) para tener una idea de la cobertura con sólo 10 bandas.