Difficulté à adapter les données linéaires morcelées dans R

Question

J'ai les données suivantes (coût d'un produit vs temps) qui ressemblent à ceci :

annum <- c(1903, 1904, 1905, 1906, 1907, 1908, 1909, 1910, 1911, 1912, 1913, 
    1914, 1915, 1916, 1917, 1918, 1919)
cost <- c(0.0000,  18.6140,  92.1278, 101.9393, 112.0808, 122.5521, 
    133.3532, 144.4843, 244.5052, 275.6068, 295.2592, 317.3145, 
    339.6527, 362.3537, 377.7775, 402.8443, 437.5539)

mydata <- as.data.frame(cbind(annum, cost))

g <- ggplot(mydata, aes(x = annum, y = cost))
g <- g + geom_point()
g <- g + scale_y_continuous(labels=scales::dollar_format())
g

Voici le graphique résultant de ces données en utilisant ce code Le graphique montre quelque chose qui ressemble à une pièce linéaire pour moi ; il y a un saut de 1904 à 1905 ; puis une ligne claire de 1905 à 1910 ; puis un saut ; et ensuite une autre ligne de 1911 à la fin. (Le premier point (1903, 0) est fictif.)

J'ai essayé d'utiliser le package segmenté pour modéliser cela, mais au lieu de choisir quelque chose comme 1904.5 et 1910.5 comme points de rupture, il trouve deux points entre 1911 et 1912.

J'ai essayé d'autres techniques (par ex. "brute force" de "The R Book," et ajustement direct), mais je ne comprends clairement pas cela autant que j'en ai besoin. Toute aide serait très appréciée.

Idéalement, j'aimerais obtenir une équation pour chaque segment et un seul graphique montrant l'ajustement en pièces et un intervalle de confiance pour l'ajustement.

Answer 1

On peut utiliser le package strucchange pour cela. Voici une version de code simplifiée :

library("strucchange")

startyear <- startyear
cost <- c(0.0000,  18.6140,  92.1278, 101.9393, 112.0808, 122.5521, 
          133.3532, 144.4843, 244.5052, 275.6068, 295.2592, 317.3145, 
          339.6527, 362.3537, 377.7775, 402.8443, 437.5539)

ts <- ts(cost, start=1903)
plot(ts)

## pour les petits ensembles de données, vous pouvez envisager de réduire la longueur du segment
bp <- breakpoints(ts ~ time(ts), data=ts, h = 5)

## Sélection des breakpoints par le BIC
plot(bp)
breakdates(bp)
fm1 <- lm(ts ~ time(ts) * breakfactor(bp), data=ts)
coef(fm1)

plot(ts, type="p")
lines(ts(fitted(fm1),  start = startyear),  col = 4)
lines(bp)
confint(bp)

lines(confint(bp))

Vous pouvez trouver plus d'informations dans la vignette du package ou l'une des publications associées, par exemple https://doi.org/10.18637/jss.v007.i02. Il est ainsi possible de réaliser des tests de signification, d'estimer des intervalles de confiance ou d'inclure des covariables.

Une longueur de segment de 2 n'est pas possible, car la variance résiduelle ne peut pas être estimée. De même, les intervalles de confiance ne peuvent être estimés que si les segments sont suffisamment longs. Par conséquent, un seul point de rupture est montré ci-dessous, alors que l'excellente réponse de @Rui Barradas omet les intervalles de confiance mais montre deux points de rupture.

Voici un exemple sans les deux premiers points et une hypothèse supplémentaire pour estimer l'intervalle de confiance en cas de petit segment :

library("strucchange")

startyear <- 1905
cost <- c(92.1278, 101.9393, 112.0808, 122.5521, 
          133.3532, 144.4843, 244.5052, 275.6068, 295.2592, 317.3145, 
          339.6527, 362.3537, 377.7775, 402.8443, 437.5539)

ts <- ts(cost, start=startyear)
bp <- breakpoints(ts ~ time(ts), data=ts, h = 5)
fm1 <- lm(ts ~ time(ts) * breakfactor(bp), data=ts)
plot(ts, type="p")
lines(ts(fitted(fm1),  start = startyear),  col = 4)
lines(confint(bp, het.err=FALSE))

Édition :

• correction des bugs de la version originale
• ajout des coefficients et de l'intervalle de confiance
• ajout des images
• exemple avec omission des deux premières valeurs ajouté

Answer 2

Voici une autre solution avec le package strucchange mais sans créer d'abord une série temporelle.

library(strucchange)

# obtenez d'abord une taille de segment comme une fraction
# du nombre d'observations
n <- nrow(mydata)
segmts <- 3
h <- (segmts + 1)/n

# estimer maintenant les points de rupture
b <- breakpoints(cost ~ annum, h = h, breaks = (segmts - 1L), data = mydata)
bp <- mydata[b$breakpoints, "annum"]

# créer une variable de regroupement pour `ggplot`
# chaque groupe est un segment
bp <- c(bp, Inf)
mydata$grp <- findInterval(mydata$annum, bp, left.open = TRUE)

# tracer les régressions linéaires
g + geom_smooth(
  mapping = aes(group = grp),
  method = "lm",
  formula = y ~ x,
  se = FALSE
)

---

Si les premiers points de données sont supprimés, il n'y aura que deux segments mais le code ci-dessus fonctionnera toujours.

mydata <- mydata[-(1:2), ]
n <- nrow(mydata)
segmts <- 2
h <- (segmts + 1)/n
b <- breakpoints(cost ~ annum, h = h, breaks = segmts - 1L, data = mydata)
bp <- mydata[b$breakpoints, "annum"]
bp <- c(bp, Inf)
mydata$grp <- findInterval(mydata$annum, bp, left.open = TRUE)
mydata$grp <- factor(mydata$grp)

g + geom_smooth(
  mapping = aes(group = grp),
  method = "lm",
  formula = y ~ x,
  se = FALSE
)

Answer 3

Les intervalles de confiance pour les problèmes de changement de points sont un problème difficile pour les méthodes fréquentistes, telles que strucchange. Souvent, vous obtenez simplement des intervalles de confiance pour chaque segment, c'est-à-dire des ruptures brutales entre les segments plutôt que des transitions fluides.

C'est plus simple en utilisant des méthodes bayésiennes. Voici une solution utilisant le package mcp. Juste pour montrer, nous affichons à la fois l'intervalle ajusté (lignes rouges en pointillés) et l'intervalle de prédiction (lignes vertes en pointillés). Les lignes grises sont des tirages aléatoires de la distribution postérieure et les densités sur l'axe des x sont les postérieurs pour les emplacements des points de changement.

data = data.frame(
  annum = 1903:1919,
  cost = c(0.0000,  18.6140,  92.1278, 101.9393, 112.0808, 122.5521, 
          133.3532, 144.4843, 244.5052, 275.6068, 295.2592, 317.3145, 
          339.6527, 362.3537, 377.7775, 402.8443, 437.5539)
)

# Modèle avec trois pentes disjointes
model = list(
  cost ~ 1 + annum,
  ~ 1 + annum,
  ~ 1 + annum
)

library(mcp)
fit = mcp(model, data)
plot(fit, q_fit = TRUE, q_predict = TRUE)

Si vous êtes intéressé par les estimations des paramètres pour les points de changement et les segments, il suffit d'appeler summary(fit):

        name    mean  lower    upper Rhat n.eff
     annum_1   -0.11   -0.2 -6.6e-04  2.5    25
     annum_2   10.36    7.4  1.3e+01  1.0   609
     annum_3   22.74   21.2  2.4e+01  1.0   264
        cp_1 1904.50 1904.0  1.9e+03  2.5    24
        cp_2 1910.46 1910.0  1.9e+03  1.0   778
 Intercept_1  221.39   10.8  3.9e+02  1.0   948
 Intercept_2   86.77   75.0  9.8e+01  1.0  1297
 Intercept_3  236.03  221.7  2.5e+02  1.0   237
     sigma_1    5.97    3.6  8.9e+00  1.0  1709

Answer 4

Cela aide-t-il. Utiliser la méthode loess?

library(tidyverse)
ggplot(mydata, aes(x = annum, y = cost))+
  geom_point()+
  geom_smooth(method = "loess", formula = "y~x")