Dapatkah seseorang tolong beri tahu saya bagaimana membuat R memperkirakan break point dalam model linear piecewise (sebagai parameter tetap atau acak), ketika saya juga perlu memperkirakan efek acak lainnya?
Saya telah memasukkan contoh mainan di bawah ini yang cocok dengan tongkat hoki / regresi tongkat patah dengan varians kemiringan acak dan varians y-intersep acak untuk titik istirahat 4. Saya ingin memperkirakan titik istirahat alih-alih menentukannya. Ini bisa berupa efek acak (lebih disukai) atau efek tetap.
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")
}
)
Keluaran:
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
terkunci
sumber
sumber
Jawaban:
Pendekatan lain adalah untuk membungkus panggilan ke lmer dalam fungsi yang melewati breakpoint sebagai parameter, kemudian meminimalkan penyimpangan dari model yang dipasang tergantung pada breakpoint menggunakan optimisasi. Ini memaksimalkan kemungkinan log profil untuk breakpoint, dan, secara umum (yaitu, tidak hanya untuk masalah ini) jika fungsi interior ke pembungkus (lebih dalam kasus ini) menemukan perkiraan kemungkinan maksimum tergantung pada parameter yang diteruskan, keseluruhan Prosedur menemukan estimasi kemungkinan maksimum bersama untuk semua parameter.
Untuk mendapatkan interval kepercayaan untuk breakpoint, Anda bisa menggunakan kemungkinan profil . Tambahkan, misalnya,
qchisq(0.95,1)
ke penyimpangan minimum (untuk interval kepercayaan 95%) kemudian cari poinfoo(x)
yang sama dengan nilai yang dihitung:Agak asimetris, tapi tidak presisi buruk untuk masalah mainan ini. Alternatifnya adalah bootstrap prosedur estimasi, jika Anda memiliki cukup data untuk membuat bootstrap dapat diandalkan.
sumber
Solusi yang diusulkan oleh jbowman sangat bagus, hanya menambahkan beberapa pernyataan teoretis:
Mengingat diskontinuitas fungsi indikator yang digunakan, kemungkinan profil mungkin sangat tidak menentu, dengan beberapa minimum lokal, jadi pengoptimal biasa mungkin tidak berfungsi. Solusi yang biasa untuk "model ambang" seperti itu adalah dengan menggunakan pencarian grid yang lebih rumit, mengevaluasi penyimpangan pada setiap hari breakpoint / ambang batas yang mungkin direalisasikan (dan tidak pada nilai di antaranya, seperti yang dilakukan dalam kode). Lihat kode di bawah.
Dalam model non-standar ini, di mana breakpoint diperkirakan, penyimpangan biasanya tidak memiliki distribusi standar. Prosedur yang lebih rumit biasanya digunakan. Lihat referensi untuk Hansen (2000) di bawah ini.
Bootstrap tidak selalu konsisten dalam hal ini, lihat Yu (akan datang) di bawah ini.
Akhirnya, tidak jelas bagi saya mengapa Anda mengubah data dengan memusatkan kembali di sekitar Days (yaitu, bp - x, bukan hanya x). Saya melihat dua masalah:
Referensi standar untuk ini adalah:
Kode:
sumber
Anda dapat mencoba model MARS . Namun, saya tidak yakin bagaimana cara menentukan efek acak.
earth(Reaction~Days+Subject, sleepstudy)
sumber
Ini adalah makalah yang mengusulkan efek campuran MARS. Seperti yang disebutkan @lockedoff, saya tidak melihat implementasi yang sama dalam paket apa pun.
sumber