¿Cómo calcula predict.lm () el intervalo de confianza y el intervalo de predicción?

Ejecuté una regresión:

CopierDataRegression <- lm(V1~V2, data=CopierData1) 

y mi tarea era obtener una

  • Intervalo de confianza del 90% para la respuesta media dada V2=6 y
  • Intervalo de predicción del 90% cuando V2=6 .

Use el siguiente código:

 X6 <- data.frame(V2=6) predict(CopierDataRegression, X6, se.fit=TRUE, interval="confidence", level=0.90) predict(CopierDataRegression, X6, se.fit=TRUE, interval="prediction", level=0.90) 

y obtuve (87.3, 91.9) y (74.5, 104.8) que parece ser correcto ya que el PI debería ser más amplio.

La salida para ambos también incluyó se.fit = 1.39 que era lo mismo. No entiendo lo que es este error estándar. ¿No debería el error estándar ser mayor para el PI frente al IC? ¿Cómo encuentro estos dos errores estándar diferentes en R? enter image description here


Datos:

 CopierData1 <- structure(list(V1 = c(20L, 60L, 46L, 41L, 12L, 137L, 68L, 89L, 4L, 32L, 144L, 156L, 93L, 36L, 72L, 100L, 105L, 131L, 127L, 57L, 66L, 101L, 109L, 74L, 134L, 112L, 18L, 73L, 111L, 96L, 123L, 90L, 20L, 28L, 3L, 57L, 86L, 132L, 112L, 27L, 131L, 34L, 27L, 61L, 77L), V2 = c(2L, 4L, 3L, 2L, 1L, 10L, 5L, 5L, 1L, 2L, 9L, 10L, 6L, 3L, 4L, 8L, 7L, 8L, 10L, 4L, 5L, 7L, 7L, 5L, 9L, 7L, 2L, 5L, 7L, 6L, 8L, 5L, 2L, 2L, 1L, 4L, 5L, 9L, 7L, 1L, 9L, 2L, 2L, 4L, 5L)), .Names = c("V1", "V2"), class = "data.frame", row.names = c(NA, -45L)) 

Al especificar un argumento de interval y level , predict.lm puede devolver el intervalo de confianza (IC) o el intervalo de predicción (IP). Esta respuesta muestra cómo obtener CI y PI sin establecer estos argumentos. Hay dos maneras:

  • utilizar el resultado de la etapa predict.lm de predict.lm ;
  • hacer todo desde cero.

Saber cómo trabajar de ambas maneras le permite comprender a fondo el procedimiento de predicción.

Tenga en cuenta que solo cubriremos el caso type = "response" (predeterminado) para predict.lm . La discusión de type = "terms" está más allá del scope de esta respuesta.


Preparar

Recojo su código aquí para ayudar a otros lectores a copiar, pegar y ejecutar. También cambio los nombres de las variables para que tengan significados más claros. Además, newdat para incluir más de una fila, para mostrar que nuestros cálculos están “vectorizados”.

 dat <- structure(list(V1 = c(20L, 60L, 46L, 41L, 12L, 137L, 68L, 89L, 4L, 32L, 144L, 156L, 93L, 36L, 72L, 100L, 105L, 131L, 127L, 57L, 66L, 101L, 109L, 74L, 134L, 112L, 18L, 73L, 111L, 96L, 123L, 90L, 20L, 28L, 3L, 57L, 86L, 132L, 112L, 27L, 131L, 34L, 27L, 61L, 77L), V2 = c(2L, 4L, 3L, 2L, 1L, 10L, 5L, 5L, 1L, 2L, 9L, 10L, 6L, 3L, 4L, 8L, 7L, 8L, 10L, 4L, 5L, 7L, 7L, 5L, 9L, 7L, 2L, 5L, 7L, 6L, 8L, 5L, 2L, 2L, 1L, 4L, 5L, 9L, 7L, 1L, 9L, 2L, 2L, 4L, 5L)), .Names = c("V1", "V2"), class = "data.frame", row.names = c(NA, -45L)) lmObject <- lm(V1 ~ V2, data = dat) newdat <- data.frame(V2 = c(6, 7)) 

A continuación, se predict.lm el resultado de predict.lm , que se comparará con nuestros cálculos manuales más adelante.

 predict(lmObject, newdat, se.fit = TRUE, interval = "confidence", level = 0.90) #$fit # fit lwr upr #1 89.63133 87.28387 91.9788 #2 104.66658 101.95686 107.3763 # #$se.fit # 1 2 #1.396411 1.611900 # #$df #[1] 43 # #$residual.scale #[1] 8.913508 predict(lmObject, newdat, se.fit = TRUE, interval = "prediction", level = 0.90) #$fit # fit lwr upr #1 89.63133 74.46433 104.7983 #2 104.66658 89.43930 119.8939 # #$se.fit # 1 2 #1.396411 1.611900 # #$df #[1] 43 # #$residual.scale #[1] 8.913508 

Use el resultado de la etapa predict.lm de predict.lm

 ## use `se.fit = TRUE` z <- predict(lmObject, newdat, se.fit = TRUE) #$fit # 1 2 # 89.63133 104.66658 # #$se.fit # 1 2 #1.396411 1.611900 # #$df #[1] 43 # #$residual.scale #[1] 8.913508 

¿Qué es se.fit ?

z$se.fit es el error estándar de la media pronosticada z$fit , utilizada para construir CI para z$fit . También necesitamos cuantiles de distribución t con un grado de libertad z$df .

 alpha <- 0.90 ## 90% Qt <- c(-1, 1) * qt((1 - alpha) / 2, z$df, lower.tail = FALSE) #[1] -1.681071 1.681071 ## 90% confidence interval CI <- z$fit + outer(z$se.fit, Qt) colnames(CI) <- c("lwr", "upr") CI # lwr upr #1 87.28387 91.9788 #2 101.95686 107.3763 

Vemos que esto concuerda con predict.lm(, interval = "confidence") .

¿Cuál es el error estándar para PI?

PI es más ancho que CI, ya que representa la varianza residual:

 variance_of_PI = variance_of_CI + variance_of_residual 

Tenga en cuenta que esto se define punto por punto. Para una regresión lineal no ponderada (como en su ejemplo), la varianza residual es igual en todas partes (conocida como homocedasticidad ), y es z$residual.scale ^ 2 . Por lo tanto, el error estándar para PI es

 se.PI <- sqrt(z$se.fit ^ 2 + z$residual.scale ^ 2) # 1 2 #9.022228 9.058082 

y PI está construido como

 PI <- z$fit + outer(se.PI, Qt) colnames(PI) <- c("lwr", "upr") PI # lwr upr #1 74.46433 104.7983 #2 89.43930 119.8939 

Vemos que esto concuerda con predict.lm(, interval = "prediction") .

observación

Las cosas son más complicadas si tiene una regresión lineal de peso, donde la varianza residual no es igual en todas partes, por lo que z$residual.scale ^ 2 debe ponderarse. Es más fácil construir PI para valores ajustados (es decir, no establece datos newdata cuando se usa type = "prediction" en predict.lm ), porque los pesos son conocidos (debe haberlos proporcionado mediante weight argumento weight al usar lm ) . Para la predicción fuera de muestra (es decir, pasa datos newdata a predict.lm ), predict.lm espera que le predict.lm cómo debe ponderarse la varianza residual. Necesita usar el argumento pred.var o los weights en predict.lm , de lo contrario recibirá una advertencia de predict.lm quejándose de información insuficiente para construir PI. Los siguientes son citados de ?predict.lm :

  The prediction intervals are for a single observation at each case in 'newdata' (or by default, the data used for the fit) with error variance(s) 'pred.var'. This can be a multiple of 'res.var', the estimated value of sigma^2: the default is to assume that future observations have the same error variance as those used for fitting. If 'weights' is supplied, the inverse of this is used as a scale factor. For a weighted fit, if the prediction is for the original data frame, 'weights' defaults to the weights used for the model fit, with a warning since it might not be the intended result. If the fit was weighted and 'newdata' is given, the default is to assume constant prediction variance, with a warning. 

Tenga en cuenta que la construcción de CI no se ve afectada por el tipo de regresión.


Haz todo desde cero

Básicamente, queremos saber cómo obtener fit , se.fit , df y residual.scale en z .

La media pronosticada puede calcularse mediante una multiplicación matriz-vector Xp %*% b , donde Xp es la matriz de predicción lineal b es el vector de coeficiente de regresión.

 Xp <- model.matrix(delete.response(terms(lmObject)), newdat) b <- coef(lmObject) yh <- c(Xp %*% b) ## c() reshape the single-column matrix to a vector #[1] 89.63133 104.66658 

Y vemos que esto concuerda con z$fit . La varianza-covarianza para yh es Xp %*% V %*% t(Xp) , donde V es la matriz de varianza-covarianza de b que puede ser calculada por

 V <- vcov(lmObject) ## use `vcov` function in R # (Intercept) V2 # (Intercept) 7.862086 -1.1927966 # V2 -1.192797 0.2333733 

La matriz de varianza-covarianza completa de yh no es necesaria para calcular CI o PI puntuales. Solo necesitamos su diagonal principal. Entonces, en lugar de hacer diag(Xp %*% V %*% t(Xp)) , podemos hacerlo de manera más eficiente a través de

 var.fit <- rowSums((Xp %*% V) * Xp) ## point-wise variance for predicted mean # 1 2 #1.949963 2.598222 sqrt(var.fit) ## this agrees with `z$se.fit` # 1 2 #1.396411 1.611900 

El grado residual de libertad está fácilmente disponible en el modelo ajustado:

 dof <- df.residual(lmObject) #[1] 43 

Finalmente, para calcular la varianza residual, use el estimador de Pearson:

 sig2 <- c(crossprod(lmObject$residuals)) / dof # [1] 79.45063 sqrt(sig2) ## this agrees with `z$residual.scale` #[1] 8.913508 

observación

Tenga en cuenta que en caso de regresión ponderada, sig2 debe calcularse como

 sig2 <- c(crossprod(sqrt(lmObject$weights) * lmObject$residuals)) / dof 

Apéndice: una función autodidacta que imita a predict.lm

El código en "Hacer todo desde cero" se ha organizado limpiamente en una función lm_predict en este Q & A: modelo lineal con lm : cómo obtener la varianza de predicción de la sum de los valores predichos .

No sé si hay una manera rápida de extraer el error estándar para el intervalo de predicción, pero siempre puede retroceder los intervalos para el SE (aunque no es un enfoque súper elegante):

 m <- lm(V1 ~ V2, data = d) newdat <- data.frame(V2=6) tcrit <- qt(0.95, m$df.residual) a <- predict(m, newdat, interval="confidence", level=0.90) cat("CI SE", (a[1, "upr"] - a[1, "fit"]) / tcrit, "\n") b <- predict(m, newdat, interval="prediction", level=0.90) cat("PI SE", (b[1, "upr"] - b[1, "fit"]) / tcrit, "\n") 

Observe que el CI SE tiene el mismo valor de se.fit .