Intervalo de previsão de inicialização

29

Existe alguma técnica de autoinicialização disponível para calcular intervalos de previsão para previsões pontuais obtidas, por exemplo, por regressão linear ou outro método de regressão (k-vizinho mais próximo, árvores de regressão etc.)?

De alguma forma, sinto que a maneira às vezes proposta de apenas inicializar a previsão do ponto (veja, por exemplo, Intervalos de previsão para regressão kNN ) não está fornecendo um intervalo de previsão, mas um intervalo de confiança.

Um exemplo em R

# STEP 1: GENERATE DATA

set.seed(34345)

n <- 100 
x <- runif(n)
y <- 1 + 0.2*x + rnorm(n)
data <- data.frame(x, y)


# STEP 2: COMPUTE CLASSIC 95%-PREDICTION INTERVAL
fit <- lm(y ~ x)
plot(fit) # not shown but looks fine with respect to all relevant aspects

# Classic prediction interval based on standard error of forecast
predict(fit, list(x = 0.1), interval = "p")
# -0.6588168 3.093755

# Classic confidence interval based on standard error of estimation
predict(fit, list(x = 0.1), interval = "c")
# 0.893388 1.54155


# STEP 3: NOW BY BOOTSTRAP
B <- 1000
pred <- numeric(B)
for (i in 1:B) {
  boot <- sample(n, n, replace = TRUE)
  fit.b <- lm(y ~ x, data = data[boot,])
  pred[i] <- predict(fit.b, list(x = 0.1))
}
quantile(pred, c(0.025, 0.975))
# 0.8699302 1.5399179

Obviamente, o intervalo de inicialização de 95% básico corresponde ao intervalo de confiança de 95%, não ao intervalo de previsão de 95%. Então, minha pergunta: como fazê-lo corretamente?

Michael M
fonte
Pelo menos no caso dos mínimos quadrados comuns, você precisará de mais do que apenas previsões pontuais; você deseja usar o erro residual estimado para também construir intervalos de previsão.
Kodiologist 31/07
@duplo: obrigado por apontar isso. A duração correta dos intervalos de previsão clássicos depende diretamente da suposição de normalidade do termo de erro; portanto, se for muito otimista, certamente também a versão com bootstrap será se for derivada a partir daí. Gostaria de saber se existe, em geral, o método de inicialização trabalhando na regressão (não necessariamente OLS).
Michael M
1
Eu acho que \ textit {inferência conforme} pode ser o que você deseja, o que permite criar intervalos de previsão baseados em reamostragem que tenham uma cobertura de amostra finita válida e que não cubram demais. Existe um bom artigo disponível em arxiv.org/pdf/1604.04173.pdf , que pode ser lido como uma introdução ao tópico, e um pacote R disponível em github.com/ryantibs/conformal .
Simon Boge Brant 10/09

Respostas:

26

O método descrito abaixo é o descrito na Seção 6.3.3 de Davidson e Hinckley (1997), Bootstrap Methods e Their Application . Graças a Glen_b e seu comentário aqui . Dado que havia várias perguntas sobre a Validação cruzada sobre esse tópico, achei que valia a pena escrever.

O modelo de regressão linear é:

YEu=XEuβ+ϵEu

Temos dados, , que usamos para estimar o β como: beta OLSEu=1,2,...,Nβ

β^OLS=(XX)-1XY

Agora, queremos prever o que será para um novo ponto de dados, já que sabemos X para ele. Esse é o problema de previsão. Vamos chamar o novo X (que sabemos) X N + 1 e o novo Y (que gostaríamos de prever), Y N + 1 . A predição habitual (se do princípio de que o ε i ii d e não correlacionada com X ) é: Y p N + 1YXXXN+1YYN+1ϵEuX

YN+1p=XN+1β^OLS

O erro de previsão feito por esta previsão é:

eN+1p=YN+1-YN+1p

Podemos reescrever esta equação como:

YN+1=YN+1p+eN+1p

Agora, já calculamos. Então, se quisermos obrigado Y N + 1 em um intervalo, digamos, 90% do tempo, tudo o que precisamos fazer é estimar de forma consistente os 5 t h e 95 t h percentis / quantis de e p N + 1 , chamada eles e 5 , e 95 , e o intervalo de previsão será [ Y p N + 1 + e 5 , Y p NYN+1pYN+15th95theN+1pe5,e95.[YN+1p+e5,YN+1p+e95]

Como estimar os quantis / percentis de ? Bem, podemos escrever: e p N + 1eN+1p

eN+1p=YN+1-YN+1p=XN+1β+ϵN+1-XN+1β^OLS=XN+1(β-β^OLS)+ϵN+1

A estratégia será coletar amostras (de um modo de inicialização) muitas vezes de e, em seguida, calcular percentis da maneira usual. Assim, talvez que vai provar 10.000 vezes a partir de e p N + 1 , e em seguida, estimar os 5 t h e 95 t h percentis como o 500 t h e 9 , 500 t h mais pequenos da amostra.eN+1peN+1p5th95th500th9,500th

Para desenhar em , que pode inicializar erros (casos seria bom, também, mas estamos assumindo erros iid de qualquer maneira). Assim, em cada replicação de inicialização, você desenha N vezes com a substituição dos resíduos ajustados pela variância (veja o próximo parágrafo) para obter ε * i , em seguida, fazer nova Y * i = X i beta OLS + ε * i , OLS, em seguida, executados no novo conjunto de dados ( Y , X )XN+1(β-β^OLS)NϵEuYEu=XEuβ^OLS+ϵEu(Y,X)para obter o dessa replicação . Por fim, este sorteio de replicação em X N + 1 ( β - β OLS ) é X N + 1 ( β OLS - β * r )βrXN+1(β-β^OLS)XN+1(β^OLS-βr)

ϵϵN+1{e1,e2,...,eN}{s1-s¯,s2-s¯,...,sN-s¯}sEu=eEu/(1-hEu)hEuEu

YN+1XXN+1

  1. YN+1p=XN+1β^OLS
  2. {s1-s¯,s2-s¯,...,sN-s¯}sEu=eEu/(1-hEu)
  3. r=1,2,...,R
    • N{ϵ1,ϵ2,...,ϵN}
    • Y=Xβ^OLS+ϵ
    • βr=(XX)-1XY
    • er=Y-Xβr
    • s-s¯
    • ϵN+1,r
    • eN+1perp=XN+1(β^OLS-βr)+ϵN+1,r
  4. 5th95theN+1pe5,e95
  5. YN+1[YN+1p+e5,YN+1p+e95]

Aqui está o Rcódigo:

# This script gives an example of the procedure to construct a prediction interval
# for a linear regression model using a bootstrap method.  The method is the one
# described in Section 6.3.3 of Davidson and Hinckley (1997),
# _Bootstrap Methods and Their Application_.


#rm(list=ls())
set.seed(12344321)
library(MASS)
library(Hmisc)

# Generate bivariate regression data
x <- runif(n=100,min=0,max=100)
y <- 1 + x + (rexp(n=100,rate=0.25)-4)

my.reg <- lm(y~x)
summary(my.reg)

# Predict y for x=78:
y.p <- coef(my.reg)["(Intercept)"] + coef(my.reg)["x"]*78
y.p

# Create adjusted residuals
leverage <- influence(my.reg)$hat
my.s.resid <- residuals(my.reg)/sqrt(1-leverage)
my.s.resid <- my.s.resid - mean(my.s.resid)


reg <- my.reg
s <- my.s.resid

the.replication <- function(reg,s,x_Np1=0){
  # Make bootstrap residuals
  ep.star <- sample(s,size=length(reg$residuals),replace=TRUE)

  # Make bootstrap Y
  y.star <- fitted(reg)+ep.star

  # Do bootstrap regression
  x <- model.frame(reg)[,2]
  bs.reg <- lm(y.star~x)

  # Create bootstrapped adjusted residuals
  bs.lev <- influence(bs.reg)$hat
  bs.s   <- residuals(bs.reg)/sqrt(1-bs.lev)
  bs.s   <- bs.s - mean(bs.s)

  # Calculate draw on prediction error
  xb.xb <- coef(my.reg)["(Intercept)"] - coef(bs.reg)["(Intercept)"] 
  xb.xb <- xb.xb + (coef(my.reg)["x"] - coef(bs.reg)["x"])*x_Np1
  return(unname(xb.xb + sample(bs.s,size=1)))
}

# Do bootstrap with 10,000 replications
ep.draws <- replicate(n=10000,the.replication(reg=my.reg,s=my.s.resid,x_Np1=78))

# Create prediction interval
y.p+quantile(ep.draws,probs=c(0.05,0.95))

# prediction interval using normal assumption
predict(my.reg,newdata=data.frame(x=78),interval="prediction",level=0.90)


# Quick and dirty Monte Carlo to see which prediction interval is better
# That is, what are the 5th and 95th percentiles of Y_{N+1}
# 
# To do it properly, I guess we would want to do the whole procedure above
# 10,000 times and then see what percentage of the time each prediction 
# interval covered Y_{N+1}

y.np1 <- 1 + 78 + (rexp(n=10000,rate=0.25)-4)
quantile(y.np1,probs=c(0.05,0.95))
Conta
fonte
Obrigado pelas explicações úteis e detalhadas. Seguindo essas linhas, acho que uma técnica geral fora do OLS (técnicas baseadas em árvore, vizinho mais próximo etc.) não estará facilmente disponível, certo?
Michael M
1
Existe um para florestas aleatórias: stats.stackexchange.com/questions/49750/…, que parece semelhante.
Bill
Xβf(X,θ)
Como você generaliza os "resíduos ajustados à variância" - a abordagem OLS se baseia na alavancagem - existe um cálculo de alavancagem para um estimador arbitrário de f (X)?
David Waterworth