Sebagai pendahuluan untuk pertanyaan tentang model linear-campuran dalam R, dan untuk berbagi sebagai referensi bagi para pecinta statistik pemula / menengah, saya memutuskan untuk memposting sebagai "gaya Tanya Jawab" yang independen, langkah-langkah yang terlibat dalam perhitungan "manual" dari koefisien dan nilai prediksi regresi linier sederhana.
Contohnya adalah dengan dataset R-built-in mtcars
,, dan akan ditetapkan sebagai mil per galon yang dikonsumsi oleh kendaraan yang bertindak sebagai variabel independen, mundur dari bobot mobil (variabel kontinu), dan jumlah silinder sebagai faktor dengan tiga level (4, 6 atau 8) tanpa interaksi.
EDIT: Jika Anda tertarik dengan pertanyaan ini, Anda pasti akan menemukan jawaban yang terperinci dan memuaskan dalam posting ini oleh Matthew Drury di luar CV .
sumber
Jawaban:
Catatan : Saya telah memposting versi yang diperluas dari jawaban ini di situs web saya .
Yakin! Kita pergi ke lubang kelinci.
Lapisan pertama adalah
lm
, antarmuka terkena programmer R. Anda dapat melihat sumbernya hanya dengan mengetiklm
di konsol R. Mayoritas (seperti sebagian besar kode tingkat produksi) sibuk memeriksa input, mengatur atribut objek, dan melempar kesalahan; tapi garis ini menonjollm.fit
adalah fungsi R lain, Anda bisa menyebutnya sendiri. Sementara denganlm
mudah bekerja dengan formula dan bingkai data, membutuhkanlm.fit
matriks, jadi itu satu tingkat abstraksi dihapus. Memeriksa sumber untuklm.fit
, lebih banyak pekerjaan, dan baris yang sangat menarik berikutSekarang kita menuju suatu tempat.
.Call
adalah cara R memanggil ke kode C. Ada fungsi C, C_Cdqrls di sumber R di suatu tempat, dan kita perlu menemukannya. Ini dia .Melihat fungsi C, sekali lagi, kami menemukan sebagian besar batas memeriksa, pembersihan kesalahan, dan pekerjaan yang sibuk. Tetapi garis ini berbeda
Jadi sekarang kita menggunakan bahasa ketiga kita, R telah memanggil C yang memanggil fortran. Ini kode fortran .
Komentar pertama menceritakan semuanya
(Menariknya, sepertinya nama rutin ini diubah di beberapa titik, tetapi seseorang lupa memperbarui komentar). Jadi kita akhirnya pada titik di mana kita dapat melakukan beberapa aljabar linier, dan benar-benar menyelesaikan sistem persamaan. Ini adalah hal yang benar-benar bagus untuk fortran, yang menjelaskan mengapa kami melewati begitu banyak lapisan untuk sampai ke sini.
Komentar tersebut juga menjelaskan apa yang akan dilakukan kode
Jadi fortran akan memecahkan sistem dengan mencari dekomposisi.Q R
Hal pertama yang terjadi, dan yang paling penting, adalah
Ini memanggil fungsi fortran
dqrdc2
pada matriks input kamix
. Apa ini?Jadi kami akhirnya berhasil linpack . Linpack adalah perpustakaan aljabar linear fortran yang telah ada sejak tahun 70-an. Aljabar linear yang paling serius akhirnya menemukan cara untuk linpack. Dalam kasus kami, kami menggunakan fungsi dqrdc2
Di sinilah pekerjaan yang sebenarnya dilakukan. Butuh satu hari penuh yang baik bagi saya untuk mencari tahu apa yang dilakukan kode ini, levelnya serendah itu. Tetapi secara umum, kami memiliki matriksX dan kami ingin memfaktorkannya ke dalam produk mana Q adalah matriks ortogonal dan R adalah matriks segitiga atas. Ini adalah hal yang cerdas untuk dilakukan, karena begitu Anda memiliki Q dan R, Anda dapat menyelesaikan persamaan linear untuk regresiX= Q R Q R Q R
sangat mudah. Memang
jadi seluruh sistem menjadi
tapi adalah segitiga atas dan memiliki peringkat yang sama seperti X t X , sehingga selama masalah kita baik yang ditimbulkan, itu adalah peringkat penuh, dan kami mungkin juga hanya memecahkan sistem berkurangR XtX
Tapi ini hal yang luar biasa. adalah segitiga atas, jadi persamaan linear terakhir di sini adalah adil , jadi penyelesaian untuk β n adalah sepele. Anda kemudian bisa naik baris, satu per satu, dan gantikan di βR βn β s yang sudah Anda ketahui, setiap kali mendapatkan persamaan linear satu variabel sederhana untuk dipecahkan. Jadi, setelah Anda memiliki dan R , semuanya runtuh menjadi apa yang disebut substitusi mundur , yang mudah. Anda dapat membaca tentang ini secara lebih rinci di sini , di mana contoh kecil eksplisit sepenuhnya dikerjakan.Q R
constant * beta_n = constant
sumber
Perhitungan langkah demi langkah aktual dalam R dijelaskan dengan indah dalam jawaban oleh Matthew Drury di utas yang sama ini. Dalam jawaban ini saya ingin berjalan melalui proses membuktikan kepada diri sendiri bahwa hasil dalam R dengan contoh sederhana dapat dicapai mengikuti aljabar linear proyeksi ke ruang kolom, dan konsep kesalahan tegak lurus (titik produk), diilustrasikan dalam posting yang berbeda , dan dijelaskan dengan baik oleh Dr. Strang di Aljabar Linier dan Penerapannya , dan mudah diakses di sini .
lm
Identik dengan:
coef(lm(mpg ~ wt + as.factor(cyl)-1))
.y_hat <- HAT %*% mpg
yang memberikan nilai identik ke:cyl <- as.factor(cyl); OLS <- lm(mpg ~ wt + cyl); predict(OLS)
:sumber
beta = solve(t(X) %*% X, t(X) %*% y)
dalam praktiknya lebih akurat daripadasolve(t(X) %*% X) %*% t(X) %*% y
.R
untuk perhitungan penting, saya ingin menyarankan agar Anda mempertimbangkan untuk mengubahnya menjadi kontribusi untuk blog kami.