Wie funktioniert vorherzusagen.lm() berechnen Sie Konfidenzintervalle und Vorhersage-Intervall?

Bei der Angabe interval und level argument predict.lm zurückkehren können Konfidenzintervall (CI) oder Prognose-Intervall (PI). Diese Antwort zeigt, so erhalten CI und PI, ohne diese Argumente. Es gibt zwei Möglichkeiten:

• Verwendung der mittleren Phase Ergebnis aus predict.lm;
• tun, alles von Grund auf.

Wissen, wie man an beide Möglichkeiten geben Ihnen eine Gründliche verstehen der Vorhersage-Verfahren.

Beachten Sie, dass wir nur die type = "response" (Standard -) Fall für predict.lm. Diskussion type = "terms" sprengt den Rahmen dieser Antwort.

Setup

Entnehme ich deinen code hier zu helfen, andere Leser zum kopieren, einfügen und ausführen. Ich auch ändern, Variablen-Namen, so dass Sie klarere Bedeutungen. Darüber hinaus erweitere ich die newdat mehr als ein Zeilen, um zu zeigen, dass unsere Berechnungen sind "vektorisiert".

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

Folgenden sind die Ausgabe von predict.lm, verglichen zu werden mit unseren manuellen Berechnungen später.

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

Verwendung der mittleren Phase Ergebnis aus 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

Was ist se.fit?

z$se.fit ist der standard-Fehler des predicted mean z$fit, verwendet, um zu konstruieren, CI z$fit. Wir müssen auch die Quantile der t-Verteilung mit einem Freiheitsgrad 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

Sehen wir, dass dies stimmt überein mit der predict.lm(, interval = "confidence").

Was ist der standard-Fehler für PI?

PI ist breiter als CI, da es Konten für die residual-Varianz:

variance_of_PI = variance_of_CI + variance_of_residual

Beachten Sie, dass dieser definierte Punkt-Weise. Für eine nicht gewichtete lineare regression (wie in deinem Beispiel), residual-Varianz ist überall gleich (bekannt als homoskedastizität), und es ist z$residual.scale ^ 2. Damit werden die standard-Fehler für PI ist

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

und PI ist gebaut wie

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

Sehen wir, dass dies stimmt überein mit der predict.lm(, interval = "prediction").

Bemerkung

Dinge werden komplizierter, wenn Sie haben ein Gewicht der linearen regression, wo die restliche Varianz ist nicht überall gleich, so dass z$residual.scale ^ 2 sollten gewichtet werden. Es ist einfacher zu konstruieren PI für angepasste Werte (das heißt, Sie setzen nicht newdata bei der Verwendung type = "prediction" im predict.lm), da die GEWICHTE bekannt sind (Sie müssen es über weight argument bei der Verwendung von lm). Für die out-of-sample prediction " (das heißt, Sie übergeben ein newdata zu predict.lm), predict.lm erwartet, dass Sie sagen, wie residual-Varianz sollten gewichtet werden. Sie müssen entweder argument pred.var oder weights im predict.lm, sonst erhalten Sie eine Warnung von predict.lm beschweren, nicht genügend Informationen für den Bau PI. Im folgenden werden zitiert aus ?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.

Beachten Sie, dass der Bau von CI ist nicht betroffen von der Art der regression.

Tun, alles von Grund auf neu

Grundsätzlich wollen wir wissen, wie zu erhalten fit, se.fit, df und residual.scale im z.

Dem predicted mean berechnet werden kann durch eine matrix-Vektor-Multiplikation Xp %*% b, wo Xp ist der lineare Prädiktor-matrix und b ist regression coefficient vector.

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

Und wir sehen, dass dies stimmt überein mit der z$fit. Die Varianz-Kovarianz-für yh ist Xp %*% V %*% t(Xp), wo V ist die Varianz-Kovarianz-matrix der b welche berechnet werden kann durch

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

Die volle Varianz-Kovarianz-matrix der yh ist nicht erforderlich, um zu berechnen, point-wise CI oder PI. Wir brauchen nur Ihre Hauptdiagonale. Also anstatt das zu tun diag(Xp %*% V %*% t(Xp)) wir tun können, ist es effizienter, über

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 

Der Verbleibende Freiheitsgrad ist leicht zugänglich in die Einbauküche Modell:

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

Schließlich, zu berechnen, residual-Varianz, verwenden Pearson Schätzer:

sig2 <- c(crossprod(lmObject$residuals)) / dof
# [1] 79.45063

sqrt(sig2)  ## this agrees with `z$residual.scale`
#[1] 8.913508

Bemerkung

Beachten Sie, dass bei der gewichteten regression sig2 sollte berechnet werden, wie

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

Anhang: eine selbst geschriebene Funktion, die imitiert predict.lm

Den code bei "Mach alles neu" wurde sauber organisiert in einer Funktion lm_predict in diesem Q & A: lineares Modell mit - lm: wie man der Vorhersage-Varianz der Summe der prognostizierten Werte.

InformationsquelleAutor 李哲源