Estimation du point de rupture dans un modèle linéaire bâton / morceau par morceaux avec des effets aléatoires dans R [code et sortie inclus]

14

Quelqu'un peut-il me dire comment R peut estimer le point de rupture dans un modèle linéaire par morceaux (en tant que paramètre fixe ou aléatoire), alors que j'ai également besoin d'estimer d'autres effets aléatoires?

J'ai inclus un exemple de jouet ci-dessous qui correspond à une régression de bâton de hockey / bâton cassé avec des variances de pente aléatoires et une variance d'ordonnée à l'origine aléatoire pour un point de rupture de 4. Je veux estimer le point de rupture au lieu de le spécifier. Il peut s'agir d'un effet aléatoire (préférable) ou d'un effet fixe.

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

Production:

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égression du bâton cassé adaptée à chaque individu

verrouillé
la source
1
Est-il possible de faire de bp un effet aléatoire?
djhocking

Réponses:

20

Une autre approche consisterait à envelopper l'appel à lmer dans une fonction à laquelle le point d'arrêt est passé en tant que paramètre, puis à minimiser la déviance du modèle ajusté conditionnelle au point d'arrêt à l'aide d'optimiser. Cela maximise la probabilité du journal de profil pour le point d'arrêt, et, en général (c'est-à-dire pas seulement pour ce problème) si la fonction à l'intérieur de l'encapsuleur (lmer dans ce cas) trouve des estimations de probabilité maximale conditionnelles au paramètre qui lui est transmis, l'ensemble trouve les estimations de maximum de vraisemblance conjointes pour tous les paramètres.

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)

Pour obtenir un intervalle de confiance pour le point d'arrêt, vous pouvez utiliser la probabilité de profil . Ajouter, par exemple, qchisq(0.95,1)à la déviance minimale (pour un intervalle de confiance à 95%) puis rechercher des points où foo(x)est égal à la valeur calculée:

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

Assez asymétrique, mais pas de mauvaise précision pour ce problème de jouet. Une alternative serait de bootstrap la procédure d'estimation, si vous avez suffisamment de données pour rendre le bootstrap fiable.

jbowman
la source
Merci - c'était très utile. Est-ce que cette technique est appelée une procédure d'estimation en deux étapes, ou a-t-elle un nom standard auquel je pourrais me référer / rechercher?
verrouillé
C'est la probabilité maximale, ou le serait si lmer maximisait la probabilité (je pense que la valeur par défaut est en fait REML, vous devez passer un paramètre REML = FALSE à lmer pour obtenir des estimations ML). juste estimé de manière imbriquée plutôt que d'un seul coup. J'ai ajouté quelques éclaircissements au début de la réponse.
jbowman
J'ai eu des problèmes d'optimisation et des CI larges lors de l'inversion de la probabilité de profil avec mes données réelles, mais j'ai obtenu des CI d'amorçage plus étroits dans mon implémentation. Envisagiez-vous un bootstrap non paramétrique avec échantillonnage avec remplacement sur les vecteurs de données des sujets? C'est-à-dire, pour les données de l'étude de sommeil, cela impliquerait un échantillonnage avec remplacement à partir des 18 vecteurs (sujets) de 10 points de données, sans rééchantillonnage dans le vecteur de données d'un sujet.
verrouillé
Oui, j'envisageais un bootstrap non paramétrique comme vous le décrivez, mais en partie parce que je ne sais pas grand-chose sur les techniques de bootstrap avancées qui peuvent (ou peuvent ne pas) être applicables. Les IC et le bootstrap basés sur la vraisemblance du profil sont tous deux asymptotiquement précis, mais il se pourrait bien que le bootstrap soit nettement meilleur pour votre échantillon.
jbowman
5

La solution proposée par jbowman est très bonne, ajoutant juste quelques remarques théoriques:

  • Étant donné la discontinuité de la fonction d'indicateur utilisée, la probabilité de profil peut être très erratique, avec plusieurs minima locaux, de sorte que les optimiseurs habituels peuvent ne pas fonctionner. La solution habituelle pour de tels "modèles de seuil" consiste à utiliser à la place la recherche de grille la plus lourde, en évaluant la déviance à chaque jour de seuil / seuil de réalisation possible (et non à des valeurs intermédiaires, comme dans le code). Voir le code en bas.

  • Dans ce modèle non standard, où le point d'arrêt est estimé, la déviance n'a généralement pas la distribution standard. Des procédures plus compliquées sont généralement utilisées. Voir la référence à Hansen (2000) ci-dessous.

  • Le bootstrap n'est pas toujours cohérent à cet égard, voir Yu (à paraître) ci-dessous.

  • Enfin, il n'est pas clair pour moi pourquoi vous transformez les données en recentrant autour des jours (c'est-à-dire bp - x au lieu de seulement x). Je vois deux problèmes:

    1. Avec cette procédure, vous créez des jours artificiels tels que 6,1 jours, 4,1 etc. Je ne sais pas comment interpréter le résultat de 6,07 par exemple, puisque vous n'avez observé que des valeurs pour les jours 6 et 7? (dans un modèle de point d'arrêt standard, toute valeur du seuil entre 6 et 7 devrait vous donner le même coefficient / écart)
    2. b1 et b2 ont la signification opposée, car pour b1 les jours diminuent, alors qu'ils augmentent pour b2? Donc, le test informel sans point d'arrêt est b1! = - b2

Les références standard pour cela sont:

  • Standard OLS: Hansen (2000) Sample Splitting and Threshold Estimation, Econometrica, Vol. 68, n ° 3. (mai 2000), pp. 575-603.
  • Modèles plus exotiques: 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 (à paraître) 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)]
Matifou
la source
0

Vous pouvez essayer un modèle MARS . Cependant, je ne sais pas comment spécifier des effets aléatoires. earth(Reaction~Days+Subject, sleepstudy)

Zach
la source
1
Merci - J'ai parcouru la documentation du package, mais il ne semble pas prendre en charge les effets aléatoires.
verrouillé
0

C'est un papier qui propose un MARS à effets mixtes. Comme @lockedoff l'a mentionné, je ne vois aucune implémentation de la même chose dans aucun package.

KarthikS
la source