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

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 ( $y$$x$$k$$(a_{1}, a_{2}, \ldots, a_{k})$

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 :$k_1$$9.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 $21.7$$-1.1$$8.5$$0.27$

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

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 = $x=15$$-1.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 × $x=9.6$$21.7057 -1.1003 \times 9.6 = 11.1428$$(9.6, 11.428)$$0.2757$$11.1428 - 0.2757 \times 9.6 = 8.496$$\beta_0 - \beta_2 k_1 = 21.7057 - 1.3760 \times 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 ( β ) $a^\top\hat\beta$$a$$a^\top\text{Var}(\hat\beta)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.

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.