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
mixed-model
lme4-nlme
change-point
piecewise-linear
verrouillé
la source
la source
Réponses:
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.
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: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.
la source
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:
Les références standard pour cela sont:
Code:
la source
Vous pouvez essayer un modèle MARS . Cependant, je ne sais pas comment spécifier des effets aléatoires.
earth(Reaction~Days+Subject, sleepstudy)
la source
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.
la source