Schätzung des Bruchpunktes in einem gebrochenen Stab / stückweise linearen Modell mit zufälligen Effekten in R [Code und Ausgabe enthalten]

Kann mir bitte jemand sagen, wie R den Bruchpunkt in einem stückweisen linearen Modell (als fester oder zufälliger Parameter) abschätzen soll, wenn ich auch andere zufällige Effekte abschätzen muss?

Im Folgenden ist ein Spielzeugbeispiel aufgeführt, das eine Hockeyschläger- / gebrochener-Schläger-Regression mit zufälligen Steigungsabweichungen und einer zufälligen y-Achsenabschnittsabweichung für einen Bruchpunkt von 4 berücksichtigt. Ich möchte den Bruchpunkt schätzen, anstatt ihn anzugeben. Dies kann ein zufälliger Effekt (bevorzugt) oder ein fester Effekt sein.

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")
        }
    )

Ausgabe:

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

Ein anderer Ansatz wäre, den Aufruf an lmer in eine Funktion zu verpacken, die den Haltepunkt als Parameter übergibt, und dann die Abweichung des angepassten Modells zu minimieren, die vom Haltepunkt abhängig ist, indem die Option "Optimize" verwendet wird. Dies maximiert die Wahrscheinlichkeit Profilprotokolls für den Haltepunkt und im Allgemeinen (dh nicht nur für dieses Problem), wenn die Funktion im Inneren des Wrappers (in diesem Fall weniger) maximale Wahrscheinlichkeitsschätzungen findet, die von dem an ihn übergebenen Parameter abhängig sind, das Ganze Mit procedure werden die gemeinsamen Maximum-Likelihood-Schätzungen für alle Parameter ermittelt.

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)

Um ein Konfidenzintervall für den Haltepunkt zu erhalten, können Sie die Profilwahrscheinlichkeit verwenden . Hinzufügen, zBqchisq(0.95,1) die minimale Abweichung (für ein 95% -Konfidenzintervall) und suche dann nach Punkten, bei denen foo(x)der berechnete Wert gleich ist:

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

Etwas asymmetrisch, aber keine schlechte Präzision für dieses Spielzeugproblem. Eine Alternative wäre das Bootstrap-Verfahren, wenn Sie über genügend Daten verfügen, um den Bootstrap zuverlässig zu machen.

Die von jbowman vorgeschlagene Lösung ist sehr gut und fügt nur einige theoretische Anmerkungen hinzu:

• Angesichts der Diskontinuität der verwendeten Indikatorfunktion kann die Profilwahrscheinlichkeit mit mehreren lokalen Minima sehr unregelmäßig sein, sodass die üblichen Optimierer möglicherweise nicht funktionieren. Die übliche Lösung für solche "Schwellenwertmodelle" besteht darin, stattdessen die umständlichere Rastersuche zu verwenden und die Abweichung an jedem möglichen realisierten Haltepunkt / Schwellentag (und nicht bei Werten dazwischen, wie im Code angegeben) zu bewerten. Siehe Code unten.

• In diesem Nicht-Standardmodell, in dem der Haltepunkt geschätzt wird, weist die Abweichung normalerweise nicht die Standardverteilung auf. In der Regel werden kompliziertere Verfahren angewendet. Siehe den Verweis auf Hansen (2000) weiter unten.

• Der Bootstrap ist in dieser Hinsicht auch nicht immer konsistent, siehe Yu (in Kürze) weiter unten.

• Schließlich ist mir nicht klar, warum Sie die Daten transformieren, indem Sie sie um die Tage herum neu zentrieren (dh bp - x statt nur x). Ich sehe zwei Probleme:

  1. Mit dieser Prozedur erstellen Sie künstliche Tage wie 6.1 Tage, 4.1 usw. Ich bin nicht sicher, wie ich das Ergebnis von 6.07 interpretieren soll, da Sie nur Werte für Tag 6 und Tag 7 beobachtet haben. (In einem Standard-Breakpoint-Modell sollte jeder Wert des Schwellenwerts zwischen 6 und 7 dieselbe Abweichung ergeben.)
  2. b1 und b2 haben die entgegengesetzte Bedeutung, da für b1 die Tage abnehmen, während sie für b2 zunehmen? Der informelle Test ohne Haltepunkt lautet also b1! = - b2

Standardreferenzen hierfür sind:

• Standard OLS: Hansen (2000) Sample Splitting and Threshold Estimation, Econometrica, Vol. 68, Nr. 3. (Mai 2000), S. 575-603.
• Weitere exotische Modelle: Lee, Seo, Shin (2011) Testen auf Schwelleneffekte in Regressionsmodellen, Journal of the American Statistical Association (Theorie und Methoden) (2011), 106, 220-231
• Ping Yu (in Vorbereitung) "The Bootstrap in Threshold Regression", Econometric Theory.

Code:

# 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)]

Sie könnten ein MARS- Modell ausprobieren . Ich bin mir jedoch nicht sicher, wie ich zufällige Effekte festlegen soll. earth(Reaction~Days+Subject, sleepstudy)

Dies ist ein Papier, das einen gemischten Effekt MARS vorschlägt. Wie @lockedoff bereits erwähnte, sehe ich in keinem Paket Implementierungen derselben.