Comment prédire.lm() calculer l'intervalle de confiance et l'intervalle de prédiction?

J'ai couru une régression:

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

et ma tâche était d'obtenir un

  • 90% intervalle de confiance pour la moyenne de la réponse donnée V2=6 et
  • 90% intervalle de prédiction à quand V2=6.

J'ai utilisé le code suivant: 

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)

et j'ai eu (87.3, 91.9) et (74.5, 104.8) qui semble être correcte, car le PI doit être plus large.

La sortie pour les deux également inclus se.fit = 1.39 qui était le même. Je ne comprends pas ce que cette erreur standard est. Ne pas le standard d'erreur plus importante pour la PI vs la CI? Comment puis-je trouver ces deux erreurs-types dans la R?
Comment prédire.lm() calculer l'intervalle de confiance et l'intervalle de prédiction?

De données:

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))
En regardant ?predict.lm, il dit: "se.fit: erreur standard de la valeur prédite signifie". "Prédit signifie" en fait on dirait qu'elle ne s'applique qu'à l'intervalle de confiance. Si vous ne voulez pas voir, il suffit de mettre se.fit = FALSE.
Je vous remercie. Je suppose que ce que je demande, c'est, comment puis-je calculer les deux std erreurs dans l'image? Afin que je puisse vérifier le calcul et savoir comment ils sont dérivés.

OriginalL'auteur Mitty | 2016-06-29

  1. 32

    Lors de la spécification de interval et level argument, predict.lm pouvez retourner intervalle de confiance (IC) ou intervalle de prédiction (IP). Cette réponse montre comment obtenir des CI et PI sans la définition de ces arguments. Il y a deux façons:

    • utilisation du moyen-résultat de l'étape de predict.lm;
    • faire tout à partir de zéro.

    De savoir comment travailler avec les deux façons de vous donner une complète compréhension de la prédiction de la procédure.

    Noter que nous ne prendrons en charge la type = "response" (par défaut) pour predict.lm. Discussion de type = "terms" est au-delà de la portée de cette réponse.

    Installation

    Je rassemble ici votre code d'aider les autres lecteurs de copier, de coller et de les exécuter. J'ai aussi changer les noms de variables, de sorte qu'ils ont plus claire de la signification. En outre, je développez le newdat inclure plusieurs lignes, pour montrer que nos calculs sont "vectorisé".

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

    Voici la sortie de predict.lm, à comparer avec notre guide des calculs plus tard.

    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

    Utilisation du moyen-résultat de l'étape 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'est-ce que se.fit?

    z$se.fit est l'erreur standard de la prédite moyenne z$fit, utilisé pour construire de CI pour z$fit. Nous avons aussi besoin de quantiles de la distribution en t avec un degré de liberté 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

    Nous voyons que ceci est en accord avec predict.lm(, interval = "confidence").

    Qu'est-ce que l'erreur standard pour les PI?

    PI est plus large que celui CI, comme il représente la variance résiduelle:

    variance_of_PI = variance_of_CI + variance_of_residual

    Noter que cela est défini au point sage. Pour un non-pondérée de régression linéaire (comme dans votre exemple), la variance résiduelle est égale partout (connu sous le nom homoscedasticity), et il est z$residual.scale ^ 2. Ainsi, la norme de l'erreur de PI est

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

    et PI est construit comme

    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

    Nous voyons que ceci est en accord avec predict.lm(, interval = "prediction").

    remarque

    Les choses sont plus compliquées si vous avez un poids de régression linéaire, où la variance résiduelle n'est pas égale partout, afin que z$residual.scale ^ 2 doit être pondérée. Il est plus facile de construire des PI pour les valeurs ajustées (qui est, vous ne définissez pas newdata lors de l'utilisation de type = "prediction" dans predict.lm), car les poids sont connus (vous devez avoir fourni via weight argument lors de l'utilisation de lm). Pour des échantillons de prédiction (qui est, vous passez un newdata à predict.lm), predict.lm attend de vous le dire combien de variance résiduelle doit être pondérée. Vous avez besoin d'utiliser l'argument pred.var ou weights dans predict.lm, sinon vous obtiendrez un message d'avertissement de predict.lm se plaindre de l'insuffisance des renseignements pour la construction de PI. Les éléments suivants sont indiqués à partir du ?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.

    Noter que la construction de l'IC n'est pas affectée par le type de régression.

    Faire tout à partir de zéro

    Fondamentalement, nous voulons savoir comment faire pour obtenir fit, se.fit, df et residual.scale dans z.

    La prédite moyenne peut être calculée par une matrice-vecteur multiplication Xp %*% b, où Xp est le prédicteur linéaire de la matrice et b est le coefficient de régression de vecteur.

    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

    Et nous voyons que ceci est en accord avec z$fit. La variance-covariance pour yh est Xp %*% V %*% t(Xp), où V est la matrice de variance-covariance de b qui peut être calculée par

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

    La pleine matrice de variance-covariance de yh n'est pas nécessaire de calculer le point de sages-CI ou PI. Nous avons seulement besoin de sa diagonale principale. Donc au lieu de faire diag(Xp %*% V %*% t(Xp)), on peut le faire plus efficacement grâce à

    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

    Résiduels degré de liberté est facilement disponible dans le modèle ajusté:

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

    Enfin, pour calculer la variance résiduelle, l'utilisation de Pearson estimateur:

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

    remarque

    Noter que dans le cas de la régression pondérée, sig2 doit être calculé comme

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

    Annexe: une auto-écrite fonction qui imite predict.lm

    Le code dans "tout Faire à partir de zéro" a été soigneusement organisée dans une fonction lm_predict dans ce Q & A: modèle linéaire avec lm: comment obtenir la prédiction de la variance de la somme des valeurs prédites.

    OriginalL'auteur 李哲源

  2. 2

    Je ne sais pas si il existe un moyen rapide pour extraire l'erreur standard de l'intervalle de prédiction, mais vous pouvez toujours backsolve les intervalles pour le SE (même si c'est pas super élégant à l'approche):

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

    Avis que la CI SE, la même valeur à partir de se.fit.

    Cela a fonctionné. Je backsolved pour SE à l'aide 89.63 + - t(0.95,43)xSE = limite Inférieure où la limite Inférieure a été 87.28 pour l'IC et 74.46 pour le PI. Le SE-CI était de 1,39 et SE PI était 9.02. Si la SE pour la prédiction de l'intervalle EST plus grand que l'intervalle de confiance. Mais je ne comprends toujours pas pourquoi la sortie de la R pour l'intervalle de prévision des listes de la se.fit = 1.39. Pourquoi n'est-il pas liste 9? Merci!!!

    OriginalL'auteur MAB