Ubah analisis titik menggunakan R's nls ()

Saya mencoba menerapkan analisis "titik perubahan", atau menggunakan regresi multifase nls()dalam R.

Ini beberapa data palsu yang saya buat . Rumus yang ingin saya gunakan agar sesuai dengan data adalah:

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

Apa yang seharusnya dilakukan adalah menyesuaikan data hingga titik tertentu dengan intersep dan kemiringan tertentu ( dan β 1 ), kemudian, setelah nilai x tertentu ( δ ), tambahkan kemiringan tersebut dengan β 2 . Itulah inti dari semuanya. Sebelum titik δ , itu akan sama dengan 0, dan $\beta_0$$\beta_1$$\delta$$\beta_2$$\delta$ akan di-nolkan.$\beta_2$

Jadi, inilah fungsi saya untuk melakukan ini:

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

Dan saya mencoba menyesuaikan model dengan cara ini

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

Saya memilih parameter awal tersebut, karena saya tahu itu adalah parameter awal, karena saya membuat data.

Namun, saya mendapatkan kesalahan ini:

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

Apakah saya baru saja membuat data yang tidak menguntungkan? Saya mencoba memasang ini pada data nyata terlebih dahulu, dan mendapatkan kesalahan yang sama, dan saya baru saja mengetahui bahwa parameter awal saya tidak cukup baik.

(Awalnya saya pikir itu bisa menjadi masalah yang dihasilkan dari fakta yang maxtidak di-vektor-kan, tapi itu tidak benar. Itu memang menyusahkan untuk bekerja dengan changePoint, karenanya modifikasi berikut:

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

Posting milis R-bantuan ini menjelaskan satu cara di mana kesalahan ini dapat terjadi: rhs dari rumus ini terlalu tinggi, sehingga mengubah dua parameter secara bersamaan memberikan kecocokan yang sama dengan data. Saya tidak bisa melihat bagaimana itu berlaku untuk model Anda, tetapi mungkin itu benar.

Bagaimanapun, Anda dapat menulis fungsi tujuan Anda sendiri dan menguranginya. Fungsi berikut memberikan kesalahan kuadrat untuk titik data (x, y) dan nilai tertentu dari parameter (struktur argumen aneh fungsi adalah untuk menjelaskan cara optimkerjanya):

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

Lalu kita katakan:

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

Dan lihat:

$par
[1] 54.53436800 -0.09283594  2.07356459 48.00000006

Perhatikan bahwa untuk data palsu saya ( x <- 40:60; data <- changePoint(x, 50, 0, 2, 48) + rnorm(21, 0, 0.5)) ada banyak maxima lokal tergantung pada nilai parameter awal yang Anda berikan. Saya kira jika Anda ingin menganggap ini serius, Anda akan memanggil pengoptimal berkali-kali dengan parameter awal acak dan memeriksa distribusi hasil.

Hanya ingin menambahkan bahwa Anda dapat melakukan ini dengan banyak paket lain. Jika Anda ingin mendapatkan perkiraan ketidakpastian di sekitar titik perubahan (sesuatu yang tidak dapat dilakukan), coba mcppaket.

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

Mari kita plot dengan interval prediksi (garis hijau). Kepadatan biru adalah distribusi posterior untuk lokasi titik perubahan:

# Plot it
plot(fit, q_predict = T)

Anda dapat memeriksa parameter individual secara lebih rinci menggunakan plot_pars(fit)dan summary(fit).