In che modo predict.lm() calcola l'intervallo di confidenza e l'intervallo di previsione?

Question

ho corse una regressione:

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

e il mio compito è stato quello di ottenere una

• 90% intervallo di di confidenza per la media risposta data V2=6 e
• 90% intervallo di previsione quando V2=6.

ho usato il seguente codice:

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) 

e ho ottenuto (87.3, 91.9) e (74.5, 104.8) che sembra essere corretta in quanto il PI dovrebbe essere più ampio.

L'uscita per entrambi comprendeva anche se.fit = 1.39 che era lo stesso. Non capisco cosa sia questo errore standard. L'errore standard non dovrebbe essere maggiore per il PI rispetto al CI? Come trovo questi due diversi errori standard in R?

---

dati:

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

Una risposta breve

## no need to specify "interval"; even no need to specify "level" 
z <- predict(CopierDataRegression, X6, se.fit=TRUE) 

L'errore standard utilizzato per la C'è:

se.CI <- z$se.fit 
# [1] 1.396411 

L'errore standard utilizzato per la PI è:

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

Per calcolare la fiducia/intervallo di predizione al 90% del livello di significatività, fare:

alpha <- qt((1-0.9)/2, df = z$df) 
# [1] -1.681071 
CI <- z$fit + c(alpha, -alpha) * se.CI 
# [1] 87.28387 91.97880 
PI <- z$fit + c(alpha, -alpha) * se.PI 
# [1] 74.46433 104.79833 

---

Una risposta più con dettagli matematici

Quando si forma un modello lineare, il modello adattato è rappresentato in forma matriciale:

y = X% *% beta.cappello

È possibile ottenere beta.hat da

beta.hat <- CopierDataRegression$coefficients 
# (Intercept)   V2 
# -0.5801567 15.0352480 

e la sua matrice di covarianza da

V <- vcov(CopierDataRegression) 
#    (Intercept)   V2 
# (Intercept) 7.862086 -1.1927966 
# V2   -1.192797 0.2333733 

Quando facciamo la previsione a matrice previsione Xp, abbiamo potuto prevedere media:

Xp <- model.matrix(~ V2, X6) 
pred <- as.numeric(Xp %*% beta.hat) 
# [1] 89.63133 

e previsione varianza:

se2 <- unname(rowSums((Xp %*% V) * Xp)) 
# [1] 1.949963 

Per il 90% di intervallo di confidenza -livello, facciamo

alpha <- qt((1-0.9)/2, df = CopierDataRegression$df.residual) 
# [1] -1.681071 
CI <- pred + c(alpha, -alpha) * sqrt(se2) 
# [1] 87.28387 91.97880 

Intervallo di previsione è un intervallo più ampio, come esso ulteriori conti per incertezza di rumore sigma2. La stima di Pearson sigma2 è

sigma2 <- sum(CopierDataRegression$residuals^2)/CopierDataRegression$df.residual 
# [1] 79.45063 

Pertanto, il 90% intervallo -livello previsione è:

PI <- pred + c(alpha, -alpha) * sqrt(se2 + sigma2) 
# [1] 74.46433 104.79833 

---

Alla fine, z <- predict(CopierDataRegression, X6, se.fit=TRUE) rendimenti

• z$fit: pred sopra;
• z$se.fit: sqrt(se2) sopra;
• z$df: df.residuals sopra;
• z$residual.scale: sqrt(sigma2) precedente.

Answer 2

Non so se c'è un modo rapido per estrarre l'errore standard per l'intervallo di previsione, ma si può sempre backsolve gli intervalli per la SE (anche se non è super elegante approccio):

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

si noti che il CI sE è lo stesso valore da se.fit.