Erro padrão de declives na regressão linear por partes com pontos de interrupção conhecidos

9

A situação

Eu tenho um conjunto de dados com um dependente uma variável independente . Quero ajustar uma regressão linear por partes contínua com pontos de interrupção conhecidos / fixos que ocorrem em . Os breakpoins são conhecidos sem incerteza, então não quero calculá-los. Então eu ajuste uma regressão (OLS) no formato Aqui está um exemplo emx k ( a 1 , a 2 , , a k ) y i = β 0 + β 1 x i + β 2 max ( x i - a 1 , 0 ) + β 3 max ( x i - a 2 , 0 ) + + Β k + 1 máx ( xyxk(a1,a2,,ak)

yi=β0+β1xi+β2max(xia1,0)+β3max(xia2,0)++βk+1max(xiak,0)+ϵi
R
set.seed(123)
x <- c(1:10, 13:22)
y <- numeric(20)
y[1:10] <- 20:11 + rnorm(10, 0, 1.5)
y[11:20] <- seq(11, 15, len=10) + rnorm(10, 0, 2)

Vamos supor que o ponto de interrupção ocorra em : 9,6k19.6

mod <- lm(y~x+I(pmax(x-9.6, 0)))
summary(mod)

Coefficients:
                    Estimate Std. Error t value Pr(>|t|)    
(Intercept)          21.7057     1.1726  18.511 1.06e-12 ***
x                    -1.1003     0.1788  -6.155 1.06e-05 ***
I(pmax(x - 9.6, 0))   1.3760     0.2688   5.120 8.54e-05 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1   1

A interceptação e a inclinação dos dois segmentos são: e para o primeiro e e para o segundo, respectivamente.- 1,1 8,5 0,2721.71.18.50.27

Ponto de interrupção


Questões

  1. Como calcular facilmente a interceptação e a inclinação de cada segmento? O modelo pode ser rememetrizado para fazer isso em um cálculo?
  2. Como calcular o erro padrão de cada inclinação de cada segmento?
  3. Como testar se duas inclinações adjacentes têm as mesmas inclinações (isto é, se o ponto de interrupção pode ser omitido)?
COOLSerdash
fonte

Respostas:

7
  1. Como calcular facilmente a interceptação e a inclinação de cada segmento?

A inclinação de cada segmento é calculada simplesmente adicionando todos os coeficientes até a posição atual. Portanto, a estimativa da inclinação em é ,.- 1.1003 + 1,3760 = 0,2757x=151.1003+1.3760=0.2757

A interceptação é um pouco mais difícil, mas é uma combinação linear de coeficientes (envolvendo os nós).

No seu exemplo, a segunda linha atende à primeira em ; portanto, o ponto vermelho está na primeira linha em . Como a segunda linha passa pelo ponto com a inclinação , sua interceptação é . Obviamente, você pode juntar essas etapas e simplificar até a interceptação para o segundo segmento = .21.7057 - 1.1003 × 9,6 = 11,1428 ( 9,6 , 11,428 ) 0,2757 11,1428 - 0,2757 × 9,6 = 8,496 β 0 - β 2 k 1 = 21.7057 - 1,3760 × 9,6x=9.621.70571.1003×9.6=11.1428(9.6,11.428)0.275711.14280.2757×9.6=8.496β0β2k1=21.70571.3760×9.6

O modelo pode ser reparameterizado para fazer isso em um cálculo?

Bem, sim, mas provavelmente é mais fácil, em geral, calculá-lo a partir do modelo.

2. Como calcular o erro padrão de cada inclinação de cada segmento?

Como a estimativa é uma combinação linear de coeficientes de regressão , onde consiste em 1 e 0s, a variação é . O erro padrão é a raiz quadrada dessa soma de termos de variância e covariância. um um Var ( β ) umaβ^aaVar(β^)a

Por exemplo, no seu exemplo, o erro padrão da inclinação do segundo segmento é:

Sb <- vcov(mod)[2:3,2:3]
sqrt(sum(Sb))

alternativamente, na forma de matriz:

Sb <- vcov(mod)
a <- matrix(c(0,1,1),nr=3)
sqrt(t(a) %*% Sb %*% a)

3. Como testar se duas inclinações adjacentes têm as mesmas inclinações (isto é, se o ponto de interrupção pode ser omitido)?

Isso é testado observando o coeficiente na tabela desse segmento. Veja esta linha:

I(pmax(x - 9.6, 0))   1.3760     0.2688   5.120 8.54e-05 ***

Essa é a mudança na inclinação em 9,6. Se essa alteração for diferente de 0, as duas inclinações não serão as mesmas. Portanto, o valor p para um teste em que o segundo segmento tem a mesma inclinação que o primeiro está logo no final dessa linha.

Glen_b -Reinstate Monica
fonte
(+1) Obrigado Glen pela sua resposta. Apenas uma pequena pergunta sobre o item 2: No meu exemplo, eu precisaria da matriz de variância-covariância de xe I(pmax(x-9.6,0)), isso está correto?
COOLSerdash
Não. Editei um exemplo explícito com base no seu exemplo. Se você quiser mais detalhes, pergunte.
Glen_b -Reinstar Monica
Muito obrigado pela edição, que esclarece um pouco para mim. Entendo isso corretamente: o erro padrão é o mesmo para cada inclinação?
COOLSerdash
11
Não. O procedimento é o mesmo, mas o valor não é. O erro padrão da inclinação do primeiro segmento está na sua tabela de regressão (0,1788). O erro padrão da inclinação do segundo segmento é 0,1160. Se tivéssemos um terceiro segmento, ele envolveria mais termos de variância-covariância em sua soma (antes que a raiz quadrada seja obtida).
Glen_b -Reinstala Monica
6

Minha abordagem ingênua, que responde à pergunta 1:

mod2 <- lm(y~I((x<9.6)*x)+as.numeric((x<9.6))+
             I((x>=9.6)*x)+as.numeric((x>=9.6))-1)
summary(mod2)

#                        Estimate Std. Error t value Pr(>|t|)    
# I((x < 9.6) * x)        -1.1040     0.2328  -4.743 0.000221 ***
# as.numeric((x < 9.6))   21.7188     1.3099  16.580 1.69e-11 ***
# I((x >= 9.6) * x)        0.2731     0.1560   1.751 0.099144 .  
# as.numeric((x >= 9.6))   8.5442     2.6790   3.189 0.005704 ** 

Mas não tenho certeza se as estatísticas (em particular os graus de liberdade) são feitas corretamente, se você fizer dessa maneira.

Roland
fonte
(+1) Muito obrigado pela sua resposta. Ele fornece uma maneira muito conveniente de calcular diretamente as intercepções e inclinações, obrigado!
COOLSerdash