Wie berechnet predict.lm() Konfidenzintervall und Vorhersageintervall?

Question

lief ich eine Regression:

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

und war meine Aufgabe, einen

• 90% Konfidenzintervall für die mittlere Antwort V2=6 gegeben und
Prädiktionsintervalls zu erhalten
• 90%, wenn V2=6.

habe ich den folgenden Code:

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) 

und ich bekam (87.3, 91.9) und (74.5, 104.8), die korrekt zu sein scheinen, da sollte das PI breiter sein.

Der Ausgang für beide enthalten auch se.fit = 1.39, die gleich war. Ich verstehe nicht, was dieser Standardfehler ist. Sollte nicht der Standardfehler für den PI gegenüber dem CI größer sein? Wie finde ich diese zwei verschiedenen Standardfehler in R?

---

Daten:

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

Eine kurze Antwort

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

Der Standardfehler für CI verwendet wird:

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

Die für PI verwendet Standardfehler:

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

Vertrauen berechnen/Prognoseintervall bei 90% Signifikanzniveau, do:

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 

---

Eine längere Antwort mit mathematischen Details

Wenn Sie ein lineares Modell passen, das angepasste Modell in Form einer Matrix dargestellt wird:

y = X% *% Beta.

Hut

Sie beta.hat von

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

und seine Kovarianzmatrix von

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

Wenn wir mit Prädiktionsmatrix Xp machen Vorhersage bekommen, haben wir Mittel vorhergesagt:

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

und Vorhersage Varianz:

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

Für 90% -Niveau Konfidenzintervall wir tun

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 

Vorhersageintervall ein breiteres Intervall, wie es weiter Konten für Unsicherheit von Lärm sigma2 ist. Die Schätzung von Pearson sigma2 ist

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

somit 90% -Niveau Prädiktionsintervall ist:

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

---

Am Ende kehrt z <- predict(CopierDataRegression, X6, se.fit=TRUE)

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

Answer 2

Ich weiß nicht, ob es eine schnelle Möglichkeit ist es, die Standardfehler für das Prognoseintervall zu extrahieren, aber Sie können die Intervalle für die SE immer backsolve (auch wenn es nicht super elegant-Ansatz):

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

Beachten sie, dass die CI SE ist der gleiche Wert von se.fit.