Szacowanie punktu przerwania w złamanym drążku / częściowym modelu liniowym z losowymi efektami w R [zawiera kod i dane wyjściowe]

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

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.

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)

#Wrapper for Mixed effects model with variable break point
foo <- function(bp)
{
  mod <- lmer(Reaction ~ b1(Days, bp) + b2(Days, bp) + (b1(Days, bp) + b2(Days, bp) | Subject), data = sleepstudy)
  deviance(mod)
}

search.range <- c(min(sleepstudy$Days)+0.5,max(sleepstudy$Days)-0.5)
foo.opt <- optimize(foo, interval = search.range)
bp <- foo.opt$minimum
bp
[1] 6.071932
mod <- lmer(Reaction ~ b1(Days, bp) + b2(Days, bp) + (b1(Days, bp) + b2(Days, bp) | Subject), data = sleepstudy)

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órych foo(x)jest równa obliczonej wartości:

foo.root <- function(bp, tgt)
{
  foo(bp) - tgt
}
tgt <- foo.opt$objective + qchisq(0.95,1)
lb95 <- uniroot(foo.root, lower=search.range[1], upper=bp, tgt=tgt)
ub95 <- uniroot(foo.root, lower=bp, upper=search.range[2], tgt=tgt)
lb95$root
[1] 5.754051
ub95$root
[1] 6.923529

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.

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:

  1. Dzięki tej procedurze tworzysz sztuczne dni, takie jak 6,1 dni, 4,1 itd. Nie jestem pewien, jak interpretować na przykład wynik 6.07, ponieważ obserwowałeś tylko wartości dla dnia 6 i dnia 7? (w standardowym modelu punktu przerwania każda wartość progu od 6 do 7 powinna dać tę samą wartość współczynnika / odchylenia)
  2. b1 i b2 mają przeciwne znaczenie, ponieważ dla b1 dni maleją, a rosną dla b2? Tak więc nieformalnym testem braku punktu przerwania jest b1! = - b2

Standardowe odniesienia do tego to:

• Standardowy OLS: podział próbki i oszacowanie progu Hansena (2000), Econometrica, t. 68, nr 3. (maj 2000), str. 575–603.
• Bardziej egzotyczne modele: Lee, Seo, Shin (2011) Testowanie efektów progowych w modelach regresji, Journal of the American Statistics Association (Theory and Methods) (2011), 106, 220-231
• Ping Yu (wkrótce) Pasek startowy w regresji progowej ”, teoria ekonometryczna.

Kod:

# Using grid search over existing values:
search.grid <- sort(unique(subset(sleepstudy, Days > search.range[1] &
Days<search.range[2], "Days", drop=TRUE)))

res <- unlist(lapply(as.list(search.grid), foo))

plot(search.grid, res, type="l")
bp_grid <- search.grid[which.min(res)]

Możesz wypróbować model MARS . Nie jestem jednak pewien, jak określić losowe efekty. earth(Reaction~Days+Subject, sleepstudy)

To jest praca, która proponuje mieszane efekty MARS. Jak wspomniano @lockedoff, nie widzę żadnych takich samych implementacji w żadnym pakiecie.