Como calcular o intervalo de confiança da interceptação x em uma regressão linear?

9

Como o erro padrão de uma regressão linear é geralmente dado para a variável de resposta, estou pensando em como obter intervalos de confiança na outra direção - por exemplo, para uma interceptação x. Sou capaz de visualizar o que pode ser, mas tenho certeza de que deve haver uma maneira direta de fazer isso. Abaixo está um exemplo em R de como visualizar isso:

set.seed(1)
x <- 1:10
a <- 20
b <- -2
y <- a + b*x + rnorm(length(x), mean=0, sd=1)

fit <- lm(y ~ x)
XINT <- -coef(fit)[1]/coef(fit)[2]

plot(y ~ x, xlim=c(0, XINT*1.1), ylim=c(-2,max(y)))
abline(h=0, lty=2, col=8); abline(fit, col=2)
points(XINT, 0, col=4, pch=4)
newdat <- data.frame(x=seq(-2,12,len=1000))

# CI
pred <- predict(fit, newdata=newdat, se.fit = TRUE) 
newdat$yplus <-pred$fit + 1.96*pred$se.fit 
newdat$yminus <-pred$fit - 1.96*pred$se.fit 
lines(yplus ~ x, newdat, col=2, lty=2)
lines(yminus ~ x, newdat, col=2, lty=2)

# approximate CI of XINT
lwr <- newdat$x[which.min((newdat$yminus-0)^2)]
upr <- newdat$x[which.min((newdat$yplus-0)^2)]
abline(v=c(lwr, upr), lty=3, col=4)

insira a descrição da imagem aqui

Marc na caixa
fonte
11
Você poderia confirmar este: library(boot); sims <- boot(data.frame(x, y), function(d, i) { fit <- lm(y ~ x, data = d[i,]) -coef(fit)[1]/coef(fit)[2] }, R = 1e4); points(quantile(sims$t, c(0.025, 0.975)), c(0, 0)). Para intervalos de previsão inversa, o arquivo de ajuda chemCal:::inverse.predictfornece a seguinte referência, que também pode ajudar a derivar um IC: Massart, LM, Vandenginste, BGM, Buydens, LMC, De Jong, S. Lewi, PJ, Smeyers-Verbeke, J. (1997 ) Handbook of Chemometrics and Qualimetrics: Part A, p. 200
Roland
11
O que você mostra no gráfico não é o IC para a interceptação. Você mostra os pontos em que as linhas de confiança inferior e superior das previsões cruzam o eixo.
Roland
11
Freqüentemente, na regressão linear, há um modelo que diz algo assim: para que os sejam tratados como aleatórios e os como fixos. Isso pode ser justificado dizendo que você está procurando uma distribuição condicional, dados os s. Na prática, se você coletar uma nova amostra, geralmente não são apenas os mas também os que mudam, sugerindo que em algumas circunstâncias eles também devem ser considerados aleatórios. Gostaria de saber se esta traz sobre a conveniência de
Yi=α+βxi+εiwhere ε1,εni.i.d. N(0,σ2),
YxxYx
Michael Hardy
11
@AdrienRenaud - Parece-me que sua resposta é excessivamente simplista, considerando os aspectos assimétricos que mencionei e são destacados pelo exercício de inicialização que Roland ilustrou. Se não estou pedindo muito, talvez você possa expandir a abordagem de probabilidade mencionada.
Marc na caixa

Respostas:

8

Como calcular o intervalo de confiança da interceptação x em uma regressão linear?

Pressupostos

  • Use o modelo de regressão simples .yi=α+βxi+εi
  • Erros têm distribuição normal condicional nos regressoresϵ|XN(0,σ2In)
  • Ajuste usando o mínimo quadrado comum

3 procedimentos para calcular o intervalo de confiança na interceptação x

Expansão de Taylor de primeira ordem

O seu modelo é com desvio padrão estimado e em e parâmetros e estimado covariância . Você resolveY=aX+bσaσbabσab

aX+b=0X=ba.

Então o desvio padrão em é dado por:σXX

(σXX)2=(σbb)2+(σaa)22σabab.

MIB

Veja o código de Marc na caixa em Como calcular o intervalo de confiança da interceptação x em uma regressão linear? .

CAPITANI-POLLASTRI

O CAPITANI-POLLASTRI fornece a função de distribuição cumulativa e a função de densidade para a razão de duas variáveis ​​aleatórias normais correlacionadas. Pode ser usado para calcular o intervalo de confiança da interceptação x em uma regressão linear. Este procedimento fornece resultados (quase) idênticos aos do MIB.

De fato, usando o quadrado mínimo ordinário e assumindo a normalidade dos erros, (verificado) e estão correlacionados (verificados).β^N(β,σ2(XTX)1)β^

O procedimento é o seguinte:

  • obtenha o estimador OLS para e .ab
  • obtenha a matriz variância-covariância e extraia .σa,σb,σab=ρσaσb
  • Suponha que e sigam uma distribuição Normal Correlacionada Bivariada, . Então a função de densidade e a Função de distribuição cumulativa de são dadas por CAPITANI-POLLASTRI.abN(a,b,σa,σb,ρ)xintercept=ba
  • Use a Função de distribuição cumulativa de para calcular quantis desejados e defina um intervalo de confiança.xintercept=ba

Comparação dos 3 procedimentos

Os procedimentos são comparados usando a seguinte configuração de dados:

  • x <- 1:10
  • a <- 20
  • b <- -2
  • y <- a + b * x + rnorm (comprimento (x), média = 0, sd = 1)

10000 amostras diferentes são geradas e analisadas usando os 3 métodos. O código (R) usado para gerar e analisar pode ser encontrado em: https://github.com/adrienrenaud/stackExchange/blob/master/crossValidated/q221630/answer.ipynb

  • MIB e CAPITANI-POLLASTRI fornecem resultados equivalentes.
  • A expansão de Taylor de primeira ordem difere significativamente dos dois outros métodos.
  • MIB e CAPITANI-POLLASTRI sofrem de subcobertura. Verificou-se que o 68% (95%) ci contém o valor verdadeiro 63% (92%) do tempo.
  • A expansão de primeira ordem de Taylor sofre de supercobertura. Verificou-se que o 68% (95%) ci contém o valor verdadeiro 87% (99%) do tempo.

Conclusões

A distribuição de interceptação x é assimétrica. Justifica um intervalo de confiança assimétrico. MIB e CAPITANI-POLLASTRI fornecem resultados equivalentes. Os CAPITANI-POLLASTRI têm uma boa justificativa teórica e fundamentam o MIB. O MIB e o CAPITANI-POLLASTRI sofrem de subcobertura moderada e podem ser usados ​​para definir intervalos de confiança.

Adrien Renaud
fonte
Obrigado por esta boa resposta. Este método implica que o erro padrão da interceptação x é simétrico? Os intervalos de previsão em minha figura implicam que esse não é o caso, e eu já vi referências a isso em outros lugares.
Marc na caixa
Sim, isso implica um intervalo simétrico. Se você deseja um assimétrico, pode usar uma probabilidade de perfil tratando os parâmetros do seu modelo como parâmetros incômodos. Mas é mais trabalho :)
Adrien Renaud
Você poderia explicar mais detalhadamente como obtém essa expressão para ? (σX/X)2
@fcop É uma expansão de Taylor. Dê uma olhada em en.wikipedia.org/wiki/Propagation_of_unertosty
Adrien Renaud
2

Eu recomendaria iniciar os resíduos:

library(boot)

set.seed(42)
sims <- boot(residuals(fit), function(r, i, d = data.frame(x, y), yhat = fitted(fit)) {

  d$y <- yhat + r[i]

  fitb <- lm(y ~ x, data = d)

  -coef(fitb)[1]/coef(fitb)[2]
}, R = 1e4)
lines(quantile(sims$t, c(0.025, 0.975)), c(0, 0), col = "blue")

plot resultante

O que você mostra no gráfico são os pontos em que o limite inferior / superior da faixa de confiança das previsões cruza o eixo. Não acho que esses sejam os limites de confiança da interceptação, mas talvez sejam uma aproximação aproximada.

Roland
fonte
Ótimo - isso já parece mais razoável do que o exemplo do seu comentário. Obrigado novamente.
Marc na caixa