Analiza punktu zmiany za pomocą R's nls ()

Próbuję zaimplementować analizę „punktu zmiany” lub regresję wielofazową nls()w R.

Oto kilka fałszywych danych, które stworzyłem . Formuła, której chcę użyć do dopasowania danych, to:

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

Powinno to polegać na dopasowaniu danych do określonego punktu z pewnym przecięciem i nachyleniem ( \ beta_0$\beta_0$ i $\beta_1$ ), a następnie, po określonej wartości x ( $\delta$ ), zwiększ nachylenie o \ beta_2$\beta_2$ . Właśnie o to chodzi w tym maksimum. Przed $\delta$ będzie równa 0, a $\beta_2$ zostanie wyzerowane.

Oto moja funkcja, aby to zrobić:

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

W ten sposób staram się dopasować model

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

Wybrałem te parametry początkowe, ponieważ wiem, że są to parametry początkowe, ponieważ utworzyłem dane.

Jednak pojawia się ten błąd:

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

Czy właśnie stworzyłem niefortunne dane? Najpierw próbowałem dopasować to do rzeczywistych danych i otrzymywałem ten sam błąd, i po prostu pomyślałem, że moje początkowe parametry początkowe nie były wystarczające.

(Na początku myślałem, że to może być problem wynikający z faktu, że maxnie jest wektorowy, ale to nie prawda To, czy zrobić to ból do pracy z Changepoint Czemu następującą modyfikacją:

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

Ten post na liście mailingowej R-help opisuje jeden ze sposobów, w jaki może wystąpić ten błąd: rh formuły jest nadparametryzowany, tak że zmiana dwóch parametrów w tandemie zapewnia takie samo dopasowanie do danych. Nie rozumiem, jak to jest w przypadku twojego modelu, ale może tak jest.

W każdym razie możesz napisać własną funkcję celu i ją zminimalizować. Poniższa funkcja podaje błąd kwadratowy dla punktów danych (x, y) i pewną wartość parametrów (dziwna struktura argumentów funkcji ma uwzględniać sposób optimdziałania):

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

Następnie mówimy:

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

I zobaczyć:

$par
[1] 54.53436800 -0.09283594  2.07356459 48.00000006

Zauważ, że w przypadku moich fałszywych danych ( x <- 40:60; data <- changePoint(x, 50, 0, 2, 48) + rnorm(21, 0, 0.5)) istnieje wiele lokalnych maksimów w zależności od podanych wartości początkowych parametrów. Podejrzewam, że jeśli chciałbyś wziąć to na poważnie, dzwoniłbyś do optymalizatora wiele razy z losowymi parametrami początkowymi i sprawdzałbyś rozkład wyników.

Chciałem tylko dodać, że możesz to zrobić z wieloma innymi pakietami. Jeśli chcesz uzyskać oszacowanie niepewności wokół punktu zmiany (coś, czego nls nie może zrobić), wypróbuj mcppakiet.

# Simulate the data
df = data.frame(x = 1:100)
df$y = c(rnorm(20, 50, 5), rnorm(80, 50 + 1.5*(df$x[21:100] - 20), 5))

# Fit the model
model = list(
  y ~ 1,  # Intercept
  ~ 0 + x  # Joined slope
)
library(mcp)
fit = mcp(model, df)

Narysujmy to z przedziałem prognozy (zielona linia). Gęstość niebieskiego jest rozkładem tylnym dla położenia punktu zmiany:

# Plot it
plot(fit, q_predict = T)

Możesz sprawdzić poszczególne parametry bardziej szczegółowo, używając plot_pars(fit)i summary(fit).