Error estándar de pendientes en regresión lineal por piezas con puntos de ruptura conocidos

Question

La situación

Tengo un conjunto de datos con un dependiente, $y$ y una de las variables independientes $x$. Quiero para adaptarse a una continua por tramos de regresión lineal con $k$ conocido/fija los breakpoints que ocurren en $(a_{1}, a_{2}, \ldots, a_{k})$. El breakpoins son conocidos, sin incertidumbre, por lo que no quiero para la estimación de ellos. Entonces me ajuste de una regresión (OLS) de la forma $$y_{i} = \beta_{0} + \beta_{1}x_{i} + \beta_{2}\operatorname{max}(x_{i}-a_{1},0) + \beta_{3}\operatorname{max}(x_{i}-a_{2},0) +\ldots+ \beta_{k+1}\operatorname{max}(x_{i}-a_{k},0) +\epsilon_{i}$$ Aquí hay un ejemplo 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 a suponer que el breakpoint $k_1$ se produce en $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

El intercepto y la pendiente de los dos segmentos son: $21.7$ $-1.1$ para el primero y $8.5$ $0.27$ para el segundo, respectivamente.

---

Preguntas

1. Cómo calcular fácilmente el intercepto y la pendiente de cada segmento? Puede el modelo reparemetrized hacer esto en este cálculo?
2. Cómo calcular el error estándar de cada una pendiente de cada segmento?
3. Cómo probar si dos laderas adyacentes tienen las mismas pendientes (es decir, si el breakpoint puede ser omitido)?

Answer 1

1. Cómo calcular fácilmente el intercepto y la pendiente de cada segmento?

La pendiente de cada segmento se calcula por la simple suma de todos los coeficientes de seguridad a la posición actual. Para la estimación de la pendiente en$x=15$$-1.1003 + 1.3760 = 0.2757\,$.

La intercepción es un poco más difícil, pero es una combinación lineal de los coeficientes de (la participación de los nudos).

En tu ejemplo, la segunda línea se reúne el primer a $x=9.6$, por lo que el punto rojo está en la primera línea en $21.7057 -1.1003 \times 9.6 = 11.1428$. Desde la segunda línea pasa a través del punto de $(9.6, 11.428)$ con pendiente $0.2757$, su intercepción es $11.1428 - 0.2757 \times 9.6 = 8.496$. Por supuesto, usted puede poner los pasos juntos y simplifica la derecha hasta la intersección de la segunda segmento = $\beta_0 - \beta_2 k_1 = 21.7057 - 1.3760 \times 9.6$.

Puede el modelo reparameterized hacer esto en este cálculo?

Bueno, sí, pero es probablemente más fácil, en general, sólo se calcula a partir del modelo.

2. Cómo calcular el error estándar de cada una pendiente de cada segmento?

Dado que la estimación es una combinación lineal de los coeficientes de regresión $a^\top\hat\beta$ donde $a$ se compone de 1 y 0, la varianza es $a^\top\text{Var}(\hat\beta)a$. El error estándar es la raíz cuadrada de la suma de la varianza y la covarianza de los términos.

por ejemplo, en tu ejemplo, el error estándar de la pendiente del segundo segmento es:

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

alternativamente, en forma de matriz:

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

3. Cómo probar si dos laderas adyacentes tienen las mismas pendientes (es decir, si el breakpoint puede ser omitido)?

Esto es probado por mirar el coeficiente de la tabla de ese segmento. Ver esta línea:

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

Ese es el cambio en la pendiente en 9.6. Si ese cambio es diferente de 0, las dos vertientes no son los mismos. Por lo que el valor p para una prueba de que el segundo segmento tiene la misma pendiente que la primera está justo al final de la línea.

Answer 2

Mi enfoque ingenuo, que responde a la pregunta 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 ** 
 

Pero no estoy seguro si las estadísticas (en particular grados de libertad) se hacen correctamente, si lo haces de esta manera.