Änderungspunktanalyse mit Rs nls ()

Ich versuche, eine "Änderungspunkt" -Analyse oder eine mehrphasige Regression mit nls()in R zu implementieren .

Hier sind einige gefälschte Daten, die ich gemacht habe . Die Formel, die ich verwenden möchte, um die Daten anzupassen, lautet:

$y = \beta_0 + \beta_1x + \beta_2\max(0,x-\delta)$

Dies soll dazu dienen, die Daten bis zu einem bestimmten Punkt mit einem bestimmten und einer bestimmten Steigung ( und ) und dann nach einem bestimmten x-Wert ( ) die Steigung um erhöhen . Darum geht es bei der ganzen Max-Sache. Vor dem Punkt ist es gleich 0 und wird auf Null gesetzt.$\beta_0$$\beta_1$$\delta$$\beta_2$$\delta$$\beta_2$

Also, hier ist meine Funktion, um dies zu tun:

changePoint <- function(x, b0, slope1, slope2, delta){ 
   b0 + (x*slope1) + (max(0, x-delta) * slope2)
}

Und ich versuche, das Modell so anzupassen

nls(y ~ changePoint(x, b0, slope1, slope2, delta), 
    data = data, 
    start = c(b0 = 50, slope1 = 0, slope2 = 2, delta = 48))

Ich habe diese Startparameter gewählt, weil ich weiß, dass dies die Startparameter sind, weil ich die Daten erfunden habe.

Ich erhalte jedoch diesen Fehler:

Error in nlsModel(formula, mf, start, wts) : 
  singular gradient matrix at initial parameter estimates

Habe ich gerade unglückliche Daten gemacht? Ich habe zuerst versucht, dies auf reale Daten anzupassen, und dabei den gleichen Fehler festgestellt. Ich habe nur herausgefunden, dass meine anfänglichen Startparameter nicht gut genug waren.

(Zuerst dachte ich, es könnte ein Problem sein, das sich aus der Tatsache ergibt, dass maxes nicht vektorisiert ist, aber das ist nicht wahr. Es macht die Arbeit mit changePoint schwierig, weshalb die folgende Änderung vorgenommen wurde:

changePoint <- function(x, b0, slope1, slope2, delta) { 
   b0 + (x*slope1) + (sapply(x-delta, function (t) max(0, t)) * slope2)
}

In diesem Beitrag zur R-Help-Mailingliste wird eine Möglichkeit beschrieben, wie dieser Fehler auftreten kann: Die rechte Hälfte der Formel ist überparametrisiert, sodass das Ändern von zwei Parametern gleichzeitig die gleiche Übereinstimmung mit den Daten ergibt. Ich kann nicht sehen, wie das für Ihr Modell zutrifft, aber vielleicht ist es das.

In jedem Fall können Sie Ihre eigene Zielfunktion schreiben und minimieren. Die folgende Funktion gibt den quadratischen Fehler für Datenpunkte (x, y) und einen bestimmten Wert der Parameter an (die seltsame Argumentstruktur der Funktion soll die Funktionsweise berücksichtigen optim):

sqerror <- function (par, x, y) {
  sum((y - changePoint(x, par[1], par[2], par[3], par[4]))^2)
}

Dann sagen wir:

optim(par = c(50, 0, 2, 48), fn = sqerror, x = x, y = data)

Und sehen:

$par
[1] 54.53436800 -0.09283594  2.07356459 48.00000006

Beachten Sie, dass es für meine gefälschten Daten ( x <- 40:60; data <- changePoint(x, 50, 0, 2, 48) + rnorm(21, 0, 0.5)) viele lokale Maxima gibt, abhängig von den anfänglichen Parameterwerten, die Sie angeben. Ich nehme an, wenn Sie dies ernst nehmen wollten, würden Sie den Optimierer oft mit zufälligen Anfangsparametern aufrufen und die Verteilung der Ergebnisse untersuchen.