Сегментированные точки останова пакетов являются переменными и обнаруживают стандартные ошибки на точках останова

Question

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

Загрузка и просмотр данных:

gbay<-data.frame(matrix(,nrow=46,ncol=3))
colnames(gbay)<-c("pop","cal.length","temp")

gbay$cal.length<-c(0.597, 0.834, 1.120, 1.353, 0.119, 1.301, 0.944, 3.127, 3.375, 3.171, 3.400, 3.376, 3.322, 3.785, 3.304, 3.197, 3.216,
 4.183, 2.171, 3.989, 3.187 ,4.153, 3.252, 4.960, 4.268, 4.827, 4.869, 3.932, 4.573, 4.645, 4.634, 4.713, 4.879, 4.724,
5.031, 4.746, 5.047, 5.714, 5.227, 4.701,4.280, 5.296, 4.977, 4.932, 4.374, 4.758)

gbay$temp<-c(16, 16, 16, 16, 16, 16, 16, 20, 20, 20, 20, 20, 20, 20, 20, 24, 24, 24, 24, 24, 24, 24, 24, 26, 26, 26, 26, 26, 26, 26, 28, 28, 28, 28,
28, 28, 28, 28, 28, 30, 30, 30, 30, 30, 30, 30)
gbay$pop<-gb

ggplot(gbay,aes(x=temp,y=cal.length))+geom_point()

1) Выше я прилагаю данные о скорости роста улитки по лабораторной температуре. Я хочу создать кусочную регрессию, которая будет определять оптимальную температуру роста и максимальную скорость роста с точкой останова. Я использовал пакет segmented, чтобы сделать это с некоторым успехом, однако у меня возникли некоторые трудности с моими данными в том, что пакет перемещается между двумя точками останова, одна на 27,56 (что, основываясь на необработанных данных и мой вопрос, компонент точки останова x, который я хочу) и другой в 18.59, который происходит с «ложной» оптимой, которую я не хочу, чтобы модель рассчитывала. Я пытался манипулировать psi в сегментированной модели, но это 50-50, какую точку останова я получаю при каждом запуске модели. Есть ли способ указать пакету сосредоточиться только на определенных границах x для поиска точки останова?

m.gbay<-glm(cal.length~temp,gbay,family=gaussian)
seg.gbay<-segmented(m.gbay,seg.Z = ~temp, psi=28)
xmin<-min(gbay$temp,na.rm=T)
xmax<-max(gbay$temp,na.rm=T)
predicted.gbay<-data.frame(temp=seq(xmin,xmax,length.out=100))
predicted.gbay$cal.length<-predict(seg.gbay,predicted.gbay)
predicted.gbay$pop<-"gb"

ggplot(predicted.gbay,aes(x=temp,y=cal.length))+geom_line(aes(x=temp,y=cal.length))+
  ylab("Shell Length (mm)")+xlab("Common Garden Temperature (°C)")

summary(seg.gbay)

2) Я пытаюсь извлечь компоненты x и y точки останова (psi) из эти данные. Я сделал это успешно. Однако я также надеялся, что мне удастся извлечь ошибку компонентов x и y для точки останова. Я думаю, что модель выдает стандартную ошибку для компонента x, но мне было интересно, есть ли способ в segmented или другом пакете найти ошибку в компоненте y точки останова?

breakpts<-data.frame(matrix(,nrow=1,ncol=4))
colnames(breakpts)<-c("brkptx","brkpty","x_err","y_err")

breakpts[1,1]<-seg.gbay$psi[[2]]
breakpts[1,2]<-(seg.gbay$psi[[2]]*coef(seg.gbay)[[2]])+(coef(seg.gbay)[[1]])
breakpts[1,3]<-seg.gbay$psi[[3]]

Answer 1

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

Я думаю, что пакет mcp делает то, что вы хотите. Вкратце, вы моделируете линию, за которой следует соединенный уклон:

model = list(
  cal.length ~ 1 + temp,  # Line with intercept
  ~ 0 + temp  # Joined slope
)

Теперь давайте подгоним ее. family = gaussian() неявно. Обратите внимание, что я установил много итераций и параллельную обработку, потому что местоположение точки изменения в этом случае довольно неидентифицируемо, поэтому для изучения апостериорного значения требуется много работы:

library(mcp)
fit = mcp(model, data = gbay, iter = 100000, cores = 3)

График по умолчанию показывает апостериор для предполагаемой точки изменения (синие линии). Мы можем добавить подогнанные интервалы (красный) и интервалы прогнозирования (зеленый). Серые линии - это 25 случайных выборок сзади:

plot(fit, q_fit = T, q_predict = T)

Как видите, задняя точка для точки изменения является бимодальной. Вы также можете позвонить plot_pars(fit) и summary(fit), чтобы увидеть более подробную информацию об отдельных параметрах. Если вы хотите проверить свидетельство для самого раннего режима по сравнению со вторым режимом, вы можете использовать hypothesis(fit, "cp_1 < 25"). Если у вас есть априорные знания о вероятных склонах, местах смены точек и т. Д. c., Это можно легко добавить с помощью, например, mcp(model, gbay, prior = list(cp_1 = "dunif(0, 25)")).

Подробнее о mcp, включая инструкции по установке, на веб-сайт mcp и этот препринт . Отказ от ответственности: я автор mcp.