Erreur standard des pentes en régression linéaire par morceaux avec des points de rupture connus

La situation

J'ai un ensemble de données avec un dépendant et une variable indépendante . Je veux adapter une régression linéaire continue par morceaux avec points d'arrêt connus / fixes se produisant à . Les breakpoins sont connus sans incertitude, donc je ne veux pas les estimer. Ensuite, j'adapte une régression (OLS) de la forme Voici un exemple dans$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)

Supposons que le point d'arrêt se produit à 9,6 :$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

L'ordonnée à l'origine et la pente des deux segments sont: et pour le premier et et pour le second, respectivement.$21.7$$-1.1$$8.5$$0.27$

Des questions

1. Comment calculer facilement l'ordonnée à l'origine et la pente de chaque segment? Le modèle peut-il être re-paramétré pour le faire en un seul calcul?
2. Comment calculer l'erreur type de chaque pente de chaque segment?
3. Comment tester si deux pentes adjacentes ont les mêmes pentes (c'est-à-dire si le point d'arrêt peut être omis)?

1. Comment calculer facilement l'ordonnée à l'origine et la pente de chaque segment?

La pente de chaque segment est calculée en ajoutant simplement tous les coefficients jusqu'à la position actuelle. Ainsi , l'estimation de la pente à est .$x=15$$-1.1003 + 1.3760 = 0.2757\,$

L'interception est un peu plus difficile, mais c'est une combinaison linéaire de coefficients (impliquant les nœuds).

Dans votre exemple, la deuxième ligne rencontre la première à , donc le point rouge est sur la première ligne à . Puisque la deuxième ligne passe par le point avec une pente de , son interception est de . Bien sûr, vous pouvez assembler ces étapes et cela simplifie jusqu'à l'interception pour le deuxième segment = .$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$

Le modèle peut-il être re-paramétré pour le faire en un seul calcul?

Eh bien, oui, mais il est probablement plus facile en général de simplement le calculer à partir du modèle.

2. Comment calculer l'erreur-type de chaque pente de chaque segment?

Étant donné que l'estimation est une combinaison linéaire de coefficients de régression , où est constitué de 1 et de 0, la variance est . L'erreur standard est la racine carrée de cette somme de termes de variance et de covariance.$a^\top\hat\beta$$a$$a^\top\text{Var}(\hat\beta)a$

par exemple dans votre exemple, l'erreur standard de la pente du deuxième segment est:

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

alternativement sous forme de matrice:

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

3. Comment tester si deux pentes adjacentes ont les mêmes pentes (c'est-à-dire si le point d'arrêt peut être omis)?

Ceci est testé en regardant le coefficient dans le tableau de ce segment. Voir cette ligne:

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

C'est le changement de pente à 9,6. Si ce changement est différent de 0, les deux pentes ne sont pas identiques. Ainsi, la valeur de p pour un test que le deuxième segment a la même pente que le premier est juste à la fin de cette ligne.

Mon approche naïve, qui répond à la question 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 ** 

Mais je ne sais pas si les statistiques (en particulier les degrés de liberté) sont faites correctement, si vous le faites de cette façon.