построение моделей пороговых / кусочных / точек изменения с 95% доверительным интервалом в R - PullRequest
Новая
       98

построение моделей пороговых / кусочных / точек изменения с 95% доверительным интервалом в R

1 голос
/ 30 мая 2020

Я хотел бы построить пороговую модель с гладкими линиями 95% доверительного интервала между линейными сегментами. Вы могли подумать, что это будет просто, но я не смог найти ответа!

Мои пороги / точки останова известны, было бы здорово, если бы существовал способ визуализировать эти данные. Я пробовал сегментированный пакет, который дает следующий график:

На графике показана пороговая модель с точкой останова 5,4. Однако доверительные интервалы между линиями регрессии не являются гладкими.

Если кто-нибудь знает какой-либо способ получить плавные (т.е. без скачка между сегментами линий) линии CI между сегментированными линиями регрессии (в идеале в ggplot), это будет удивительный. Огромное спасибо.

Я включил образцы данных и код, который я пробовал ниже: 

x <- c(2.26, 1.95, 1.59, 1.81, 2.01, 1.63, 1.62, 1.19, 1.41, 1.35, 1.32, 1.52, 1.10, 1.12, 1.11, 1.14, 1.23, 1.05, 0.95, 1.30, 0.79,
0.81, 1.15, 1.10, 1.29, 0.97, 1.05, 1.05, 0.84, 0.64, 0.80, 0.81, 0.61, 0.71, 0.75, 0.30, 0.30, 0.49, 1.13, 0.55, 0.77, 0.51,
0.67, 0.43, 1.11, 0.29, 0.36, 0.57, 0.02, 0.22, 3.18, 3.79, 2.49, 2.44, 2.12, 2.45, 3.22, 3.44, 3.86, 3.53, 3.13)

y <- c(22.37, 18.93, 16.99, 15.65, 14.62, 13.79, 13.09, 12.49, 11.95, 11.48, 11.05, 10.66, 10.30,  9.96,  9.65,  9.35,  9.07,  8.81,
       8.56,  8.32,  8.09,  7.87,  7.65,  7.45,  7.25,  7.05,  6.86,  6.68,  6.50,  6.32,  6.15,  5.97,  5.80,  5.63,  5.47,  5.30,
        5.13,  4.96,  4.80,  4.63,  4.45,  4.28,  4.09,  3.90,  3.71,  3.50,  3.27,  3.01,  2.70,  2.28, 22.37, 16.99, 11.05,  8.81,
       8.56,  8.32,  7.25,  7.05,  6.50,  6.15,  5.63)

lin.mod <- lm(y ~  x)
segmented.mod <- segmented(lin.mod, seg.Z = ~x, psi=2)
plot(x, y)
plot(segmented.mod, add=TRUE, conf.level = 0.95)

, который дает следующий график (и связанные с ним скачки с 95% доверительным интервалом):

сегментированный участок

1 Ответ

2 голосов
/ 31 мая 2020

Предпосылки: Негладкость в существующих пакетах точек изменения связана с тем, что частотные пакеты работают с фиксированным значением точки изменения. Но, как и в случае со всеми предполагаемыми параметрами, это неверно, потому что действительно существует неопределенность в отношении местоположения изменения.

Решение: AFAIK, только байесовские методы могут количественно оценить это, и пакет mcp заполняет это пространство. 

library(mcp)
model = list(
  y ~ 1 + x,   # Segment 1: Intercept and slope
  ~ 0 + x  # Segment 2: Joined slope (no intercept change)
)
fit = mcp(model, data = data.frame(x, y))

График по умолчанию (plot.mcpfit() возвращает ggplot объект): 

plot(fit) + ggtitle("Default plot")

Default mcp plot

Каждая строка представляет возможную модель, которая сгенерировала данные. Апостериор для точки изменения отображается синим цветом. Вы можете добавить надежный интервал вверху, используя plot(fit, q_fit = TRUE), или построить его отдельно: 

plot(fit, lines = 0, q_fit = c(0.025, 0.975), cp_dens = FALSE) + ggtitle("Credible interval only")

mcp credible interval plot

Если ваша точка изменения действительно известно , и если вы хотите смоделировать разные остаточные масштабы для каждого сегмента (например, квази-эмулировать segmented), вы можете сделать: 

model2 = list(
  y ~ 1 + x,
  ~ 0 + x + sigma(1)  # Add intercept change in residual scale
)
fit = mcp(model2, df, prior = list(cp_1 = 1.9))  # Note: prior is a fixed value - not a distribution.
plot(fit, q_fit = TRUE, cp_dens = FALSE)

mcp emulating segmented

Обратите внимание, что CI не «прыгает» вокруг точки изменения, как в segmented. Я считаю, что это правильное поведение. Раскрытие информации: я являюсь автором mcp.

...