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

regression linear-regression

Question

Corrí una regresión:

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

y mi tarea era obtener un

• 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 IP debería ser más amplio.

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

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))

Answer 1

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

• use el resultado de la etapa predict.lm de predict.lm ;
• haz todo desde cero.

Saber cómo trabajar en ambos sentidos le brinda una comprensión profunda del 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 alcance de esta respuesta.

Preparar

Recopilo 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 el 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))

Los siguientes son los resultados de predict.lm , para compararlos 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

Utilice 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 amplio que CI, ya que explica la varianza residual:

variance_of_PI = variance_of_CI + variance_of_residual

Tenga en cuenta que esto se define puntualmente. 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 se construye 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 los valores ajustados (es decir, no establece newdata cuando usa type = "prediction" en predict.lm ), porque los pesos son conocidos (debe haberlo proporcionado a través del argumento de weight al usar lm ) . Para la predicción fuera de la muestra (es decir, si pasa datos newdata a predict.lm ), predict.lm espera que le diga cómo se debe ponderar 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. Lo siguiente se cita 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 predictora 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 calcularse mediante

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 puntualmente. 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 predict.lm que imita 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 variación de predicción de la suma de valores pronosticados .

Answer 2

No sé si hay una forma rápida de extraer el error estándar para el intervalo de predicción, pero siempre puede resolver 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 es el mismo valor de se.fit .