Pasang garis regresi yang kuat menggunakan MM-estimator di R

8

Konteks. Saya ingin mencocokkan garis regresi untuk mempelajari hubungan antara beberapa variabel respony dan beberapa kovariat terus menerus x. Karena adanya poin leverage yang buruk, saya telah memilih untuk estimator MM daripada estimasi LS.

Metodologi. Pada dasarnya, estimasi-MM adalah estimasi-M yang diinisialisasi oleh S-estimator. Oleh karena itu, dua fungsi kerugian harus dipilih. Saya telah memilih fungsi kehilangan Tukey Biweight yang banyak digunakan

ρ(u)={1[1(uk)2]3if |u|k1if |u|>k,

dengan k=1.548 di penduga-S awal (yang memberikan titik rincian sama dengan 50%), dan dengan k=2.697 pada langkah estimasi-M (untuk menjamin 70% Efisiensi Gaussian).

Saya ingin menggunakan R agar sesuai dengan garis regresi saya yang kuat.

Pertanyaan.

library(MASS)
rlm(y~x, 
    method="MM",
    k0=1.548, c=2.697,
    maxit=50)
  • Apakah kode saya konsisten dengan paragraf sebelumnya?
  • Apakah Anda akan menggunakan argumen opsional lainnya?

EDIT. Setelah diskusi saya dengan @Jason Morgan, saya menyadari bahwa kode saya sebelumnya salah. (@Jason Morgan: Terima kasih banyak untuk ini!) Namun, saya masih tidak yakin dengan lamarannya. Sebaliknya, inilah yang saya usulkan sekarang:

library(robustbase)
lmrob(y~x, 
      tuning.chi=1.548, tuning.psi=2.697)

Saya pikir itu melekat pada metodologi sekarang. Apa kamu setuju?

Terima kasih!

okram
sumber

Jawaban:

5

Secara default, dokumentasi menunjukkan bahwa rlmmenggunakan psi=psi.huberbobot. Jadi, jika Anda ingin menggunakan bisquare Tukey, Anda perlu menentukan psi=psi.bisquare. Pengaturan defaultnya adalah psi.bisquare(u, c = 4.685, deriv = 0), yang dapat Anda ubah sesuai keinginan. Misalnya, mungkin sesuatu seperti

rlm(x ~ y, method="MM", psi=psi.bisquare, maxit=50)

Anda mungkin juga ingin menyelidiki apakah Anda harus menggunakan kuadrat-terkecil ( init="lts") untuk menginisialisasi nilai awal Anda. Standarnya adalah menggunakan kuadrat terkecil.

Jason Morgan
sumber
@Janson Morgan: Anda yakin dengan apa yang Anda kemukakan? Apakah Anda punya pengalaman dengan fungsi itu? Dokumentasi saya (R 2.13.1) sebenarnya menunjukkan "Set awal koefisien dan skala akhir dipilih oleh S-estimator dengan k0 = 1,548; ini memberikan (untuk n >> p) titik kerusakan 0,5. Estimasi akhir adalah Penaksir-M dengan skala berat dan tetap Tukey yang akan mewarisi titik gangguan ini asalkan c> k0; ini berlaku untuk nilai default c yang sesuai dengan efisiensi relatif 95% pada normal. "
ocram
1
Saya telah memperkirakan model-model ini di masa lalu. Seperti yang dinyatakan dalam dokumentasi, langkah pertama dalam estimasi MM dilakukan dengan bobot Huber, yang kedua dengan bobot bisquared. Catatan saya (dari beberapa tahun yang lalu) menyatakan bahwa pada langkah-S pertama, Anda dapat menggunakan bobot bisquare alih-alih bobot Huber jika Anda menentukannya psi. Saya mungkin akan meninggalkan cdefaultnya untuk memulai (saya akan mengubah jawaban saya sesuai).
Jason Morgan
1
Saya juga menggunakan rlm, dan menggunakan fungsi psi bisquare karena properti redescending-nya. Kadang-kadang ada masalah konvergensi dengan itu, terutama dengan sampel yang lebih kecil.
jbowman