Estimación del punto de ruptura en un modelo lineal roto / a trozos con efectos aleatorios en R [código y salida incluidos]

Question

¿Puede alguien decirme cómo hacer que R estime el punto de ruptura en un modelo lineal a trozos (como parámetro fijo o aleatorio), cuando también necesito estimar otros efectos aleatorios?

He incluido un ejemplo de juguete a continuación que ajusta una regresión de palo de hockey / palo roto con varianzas aleatorias de la pendiente y una varianza aleatoria de la intersección y para un punto de ruptura de 4. Quiero estimar el punto de ruptura en lugar de especificarlo. Podría ser un efecto aleatorio (preferible) o un efecto fijo.

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

La salida:

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

Answer 1

Otro enfoque sería envolver la llamada a lmer en una función a la que se le pasa el punto de ruptura como parámetro, y luego minimizar la desviación del modelo ajustado condicionado al punto de ruptura utilizando optimize. Esto maximiza la perfil log probabilidad para el punto de ruptura, y, en general (es decir, no sólo para este problema) si la función interior al wrapper (lmer en este caso) encuentra estimaciones de máxima verosimilitud condicionadas al parámetro que se le pasa, todo el procedimiento encuentra las estimaciones conjuntas de máxima verosimilitud para todos los parámetros.

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)

Para obtener un intervalo de confianza para el punto de ruptura, puede utilizar la función perfil de probabilidad . Añade, por ejemplo, qchisq(0.95,1) a la desviación mínima (para un intervalo de confianza del 95%) y luego buscar los puntos en los que foo(x) es igual al valor calculado:

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

Algo asimétrico, pero no está mal la precisión para este problema de juguete. Una alternativa sería hacer un bootstrap del procedimiento de estimación, si se tienen suficientes datos para que el bootstrap sea fiable.

Answer 2

La solución propuesta por jbowman es muy buena, solo añado algunas observaciones teóricas:

• Dada la discontinuidad de la función indicadora utilizada, la probabilidad del perfil podría ser muy errática, con múltiples mínimos locales, por lo que los optimizadores habituales podrían no funcionar. La solución habitual para este tipo de "modelos de umbral" es utilizar en su lugar la más engorrosa búsqueda en cuadrícula, evaluando la desviación en cada uno de los posibles días de ruptura/umbral realizados (y no en los valores intermedios, como se hace en el código). Véase el código en la parte inferior.

• Dentro de este modelo no estándar, donde se estima el punto de ruptura, la desviación no suele tener la distribución estándar. Se suelen utilizar procedimientos más complicados. Véase la referencia a Hansen (2000) más abajo.

• El bootstrap tampoco es siempre consistente en este sentido, véase Yu (de próxima aparición) más adelante.

• Por último, no me queda claro por qué transformas los datos volviendo a centrarlos alrededor de los Días (es decir, bp - x en lugar de sólo x). Veo dos problemas:

  1. Con este procedimiento, se crean días artificiales como 6,1 días, 4,1, etc. No sé cómo interpretar el resultado de 6,07, por ejemplo, ya que sólo se han observado valores para el día 6 y el día 7 (en un modelo de punto de ruptura estándar, cualquier valor del umbral entre 6 y 7 debería dar el mismo coeficiente/desviación).
  2. ¿b1 y b2 tienen el significado opuesto, ya que para b1 los días disminuyen, mientras que aumentan para b2? Así que la prueba informal de que no hay punto de ruptura es b1 != - b2

Las referencias estándar para esto son:

• OLS estándar: Hansen (2000) Sample Splitting and Threshold Estimation, Econometrica, Vol. 68, No. 3. (mayo de 2000), pp. 575-603.
• Modelos más exóticos: Lee, Seo, Shin (2011) Testing for threshold effects in regression models, Journal of the American Statistical Association (Theory and Methods) (2011), 106, 220-231.
• Ping Yu (de próxima aparición) The Bootstrap in Threshold Regression", Econometric Theory.

Código:

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

Answer 3

Podrías probar con un MARS modelo. Sin embargo, no estoy seguro de cómo especificar los efectos aleatorios. earth(Reaction~Days+Subject, sleepstudy)

Answer 4

Este es un trabajo que propone un MARS de efectos mixtos. Como mencionó @lockedoff, no veo ninguna implementación del mismo en ningún paquete.