Может кто-нибудь сказать мне, как заставить R оценить точку разрыва в кусочно-линейной модели (как фиксированный или случайный параметр), когда мне также нужно оценить другие случайные эффекты?
Ниже я привел пример с игрушкой, который соответствует регрессии хоккейной клюшки / сломанной палки со случайными отклонениями наклона и случайной дисперсией y-точки пересечения для точки разрыва 4. Я хочу оценить точку разрыва вместо ее определения. Это может быть случайный эффект (предпочтительно) или фиксированный эффект.
library(lme4)
str(sleepstudy)
#Basis functions
bp = 4
b1 <- function(x, bp) ifelse(x < bp, bp - x, 0)
b2 <- function(x, bp) ifelse(x < bp, 0, x - bp)
#Mixed effects model with break point = 4
(mod <- lmer(Reaction ~ b1(Days, bp) + b2(Days, bp) + (b1(Days, bp) + b2(Days, bp) | Subject), data = sleepstudy))
#Plot with break point = 4
xyplot(
Reaction ~ Days | Subject, sleepstudy, aspect = "xy",
layout = c(6,3), type = c("g", "p", "r"),
xlab = "Days of sleep deprivation",
ylab = "Average reaction time (ms)",
panel = function(x,y) {
panel.points(x,y)
panel.lmline(x,y)
pred <- predict(lm(y ~ b1(x, bp) + b2(x, bp)), newdata = data.frame(x = 0:9))
panel.lines(0:9, pred, lwd=1, lty=2, col="red")
}
)
Выход:
Linear mixed model fit by REML
Formula: Reaction ~ b1(Days, bp) + b2(Days, bp) + (b1(Days, bp) + b2(Days, bp) | Subject)
Data: sleepstudy
AIC BIC logLik deviance REMLdev
1751 1783 -865.6 1744 1731
Random effects:
Groups Name Variance Std.Dev. Corr
Subject (Intercept) 1709.489 41.3460
b1(Days, bp) 90.238 9.4994 -0.797
b2(Days, bp) 59.348 7.7038 0.118 -0.008
Residual 563.030 23.7283
Number of obs: 180, groups: Subject, 18
Fixed effects:
Estimate Std. Error t value
(Intercept) 289.725 10.350 27.994
b1(Days, bp) -8.781 2.721 -3.227
b2(Days, bp) 11.710 2.184 5.362
Correlation of Fixed Effects:
(Intr) b1(D,b
b1(Days,bp) -0.761
b2(Days,bp) -0.054 0.181
r
mixed-model
lme4-nlme
change-point
piecewise-linear
lockedoff
источник
источник
Ответы:
Другой подход заключается в том, чтобы обернуть вызов lmer в функцию, в которой в качестве параметра передается точка останова, а затем минимизировать отклонение подобранной модели, обусловленной точкой останова, с помощью оптимизации. Это максимизирует вероятность регистрации профиля для точки останова, и, в общем (то есть, не только для этой проблемы), если внутренняя часть функции для оболочки (в данном случае lmer) находит оценки максимального правдоподобия, зависящие от передаваемого ей параметра, все Процедура находит совместные оценки максимального правдоподобия для всех параметров.
Чтобы получить доверительный интервал для точки останова, вы можете использовать профиль вероятности . Добавьте, например,
qchisq(0.95,1)
к минимальному отклонению (для доверительного интервала 95%), затем найдите точки, гдеfoo(x)
равно вычисленному значению:Несколько асимметричная, но неплохая точность для этой игрушечной задачи. Альтернативой может быть начальная процедура оценки, если у вас достаточно данных для обеспечения надежности начальной загрузки.
источник
Решение, предложенное jbowman, очень хорошее, просто добавив несколько теоретических замечаний:
Учитывая прерывистость используемой функции индикатора, вероятность профиля может быть очень ошибочной, с несколькими локальными минимумами, поэтому обычные оптимизаторы могут не работать. Обычное решение для таких «пороговых моделей» состоит в том, чтобы вместо этого использовать более громоздкий поиск в сетке, оценивая отклонение в каждой возможной реализованной точке останова / пороговом дне (а не в значениях между ними, как это сделано в коде). Смотрите код внизу.
В рамках этой нестандартной модели, где оценивается точка останова, отклонение обычно не имеет стандартного распределения. Обычно используются более сложные процедуры. См. Ссылку на Хансена (2000) ниже.
Бутстрап не всегда последовательн в этом отношении, см. Yu (готовится к печати) ниже.
Наконец, мне не ясно, почему вы преобразуете данные, перецентрируясь вокруг Дней (то есть, bp - x вместо просто x). Я вижу две проблемы:
Стандартные ссылки для этого:
Код:
источник
Вы можете попробовать модель MARS . Однако я не уверен, как указать случайные эффекты.
earth(Reaction~Days+Subject, sleepstudy)
источник
Это документ, который предлагает смешанные эффекты MARS. Как уже упоминалось @lockedoff, я не вижу ни одной реализации в одном пакете.
источник