Czy ktoś może mi powiedzieć, jak R oszacować punkt przerwania w częściowym modelu liniowym (jako parametr stały lub losowy), gdy muszę również oszacować inne efekty losowe?
Poniżej zamieściłem przykład zabawki, który pasuje do regresji kija hokejowego / łamanego kija z losowymi wariancjami nachylenia i losową wariancją przechwytywania y dla punktu złamania 4. Chcę oszacować punkt przerwania zamiast go określać. Może to być efekt losowy (najlepiej) lub efekt stały.
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")
}
)
Wynik:
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
zablokowane
źródło
źródło
Odpowiedzi:
Innym podejściem byłoby zawinięcie wywołania lmer w funkcję, której parametr przerwano jako parametr, a następnie zminimalizowanie odchylenia dopasowanego modelu w zależności od punktu przerwania przy użyciu funkcji optymalizacji. Maksymalizuje to prawdopodobieństwo dziennika profilu dla punktu przerwania i ogólnie (tj. Nie tylko w przypadku tego problemu), jeśli funkcja wnętrza opakowania (w tym przypadku lżejsza) znajdzie oszacowania maksymalnego prawdopodobieństwa zależne od przekazanego mu parametru, całość procedura wyszukuje łączne oszacowania maksymalnego prawdopodobieństwa dla wszystkich parametrów.
Aby uzyskać przedział ufności dla punktu przerwania, możesz użyć prawdopodobieństwa profilu . Dodaj np.
qchisq(0.95,1)
Do minimalnego odchylenia (dla 95% przedziału ufności), a następnie wyszukaj punkty, w którychfoo(x)
jest równa obliczonej wartości:Nieco asymetryczna, ale niezła precyzja dla tego problemu z zabawkami. Alternatywą byłoby uruchomienie procedury szacowania, jeśli masz wystarczającą ilość danych, aby uczynić bootstrap niezawodnym.
źródło
Rozwiązanie zaproponowane przez jbowman jest bardzo dobre, wystarczy dodać kilka uwag teoretycznych:
Biorąc pod uwagę nieciągłość zastosowanej funkcji wskaźnika, prawdopodobieństwo profilu może być bardzo zmienne, z wieloma lokalnymi minimami, więc zwykłe optymalizatory mogą nie działać. Zwykle rozwiązaniem dla takich „modeli progowych” jest użycie bardziej kłopotliwego wyszukiwania siatki, oceniając odchylenie w każdym możliwym zrealizowanym dniu punktu przerwania / progu (a nie wartości pośrednich, jak to zrobiono w kodzie). Zobacz kod na dole.
W tym niestandardowym modelu, w którym szacuje się punkt przerwania, odchylenie zwykle nie ma rozkładu standardowego. Zwykle stosuje się bardziej skomplikowane procedury. Zobacz odniesienie do Hansena (2000) poniżej.
Bootstrap nie zawsze jest spójny pod tym względem, patrz Yu (wkrótce) poniżej.
Wreszcie, nie jest dla mnie jasne, dlaczego transformujesz dane poprzez ponowne centrowanie wokół Dni (tj. Bp - x zamiast tylko x). Widzę dwa problemy:
Standardowe odniesienia do tego to:
Kod:
źródło
Możesz wypróbować model MARS . Nie jestem jednak pewien, jak określić losowe efekty.
earth(Reaction~Days+Subject, sleepstudy)
źródło
To jest praca, która proponuje mieszane efekty MARS. Jak wspomniano @lockedoff, nie widzę żadnych takich samych implementacji w żadnym pakiecie.
źródło