Cómo generar bandas de confianza para$\hat{Y}$

Question

Supongamos que tengo un modelo de regresión lineal. Estoy interesado en la generación de intervalos de predicción. Los valores pronosticados son fáciles de calcular, pero ¿cómo puedo calcular las desviaciones estándar para cada una de las $\hat{Y}$s?

En R, yo uso el predict.lm de la función y el uso de la interval='prediction' argumento. Yo wouldsubtract la más alta cota de la CI para el número real y dividir por 1.96 95% CI, pero me gustaría obtener directamente para estar seguro.

Traté de se.ajuste = TRUE y cosas por el estilo, pero no de trabajo.

Answer 1

Mediante el ajuste de una película objeto de obtener todos los componentes necesarios para ello. Matemáticamente se tienen estimaciones:

$\hat{\beta} = \left( \mathbf{X}^T\mathbf{X} \right) ^{-1} \left( \mathbf{X}^T y \right) $

y y de la estimación:

$$\mbox{vcov}\left(\hat{\beta} \right) = \hat{\sigma}^2 \left( \mathbf{X}^T\mathbf{X} \right) ^{-1} $$

the beta-hats are obtained by calling coef to the lm object and the variance estimate, vcov to the lm object.

Mathematically, for any $\mathbf{X}_{pred}$ observation you wish to predict the fitting $\hat{Y} = E \left[ Y | \mathbf{X} = \mathbf{X}_{pred} \right]$ then since the $\hat{Y}$ is given by: $\mathbf{X}_{pred}^T \hat{\beta}$ it is a simple mathematical manipulation to find that:

$$\mbox{var} \left( \hat{Y} \right) = \mathbf{X}_{pred}^T \mbox{vcov}\left(\hat{\beta} \right) \mathbf{X}_{pred} = \hat{\sigma}^2 \left( \mathbf{X}_{pred}^T\mathbf{X}_{pred} \right)$$

It is a simple rule of quadratic forms that the farther $\mathbf{X}_{pred}$ is from the sample mean for each covariate (in a Euclidean sense), the greater $\left( \mathbf{X}_{pred}^T\mathbf{X}_{pred} \right)$ will be and, hence, the greater the variance of $\hat{Y}$.

Simply, the variance only differs as a function of the cross product of your predicted $X$. Un ejemplo que ilustra en R ya que parecen estar interesados en los aspectos teóricos y aspectos computacionales...

x <- 1:100
y <- rnorm(100, x, 100)
plot(x, y)
f <- lm(y ~ x)
X <- model.matrix(f)
pred.se <- apply(X, 1, function(Xrow) t(Xrow) %*% vcov(f) %*% Xrow)
lines(1:100, 1:100 + 1.96*sqrt(pred.se))
lines(1:100, 1:100 - 1.96*sqrt(pred.se))
## "conf band is for uncertainty in predicted ys, should be substantially 
## tighter than observed vales

Answer 2

tenemos $$ \ hat {\ beta} \ pm t _ {\ alpha / 2, n-2} \ sqrt {\ frac {MSE} {\ sum (x_i- \ bar {x}) ^ 2}} $$ luego

   l=lm(y~x)
  MSE=mean ( (l$residuals)^2) 
      SSX=sum ( (x-mean(x))^2 )
      U= l$coefficients + qt(1-alpha/2,n-2) * sqrt(MSE/SSX)
  L= l$coefficients - qt(1-alpha/2,n-2) * sqrt(MSE/SSX)