Regresi komponen utama di r (langkah demi langkah)

Mengingat sekumpulan variabel prediktor p dan variabel respons, regresi linier berganda menggunakan metode yang dikenal sebagai kuadrat terkecil untuk meminimalkan jumlah sisa kuadrat (RSS):

RSS = Σ(y _saya – ŷ _saya ) ²

Emas:

• Σ : Simbol Yunani yang berarti jumlah
• y _i : nilai respon sebenarnya untuk observasi ^ke-i
• ŷ _i : Nilai respons yang diprediksi berdasarkan model regresi linier berganda

Namun, ketika variabel prediktor berkorelasi tinggi, multikolinearitas bisa menjadi masalah. Hal ini dapat membuat estimasi koefisien model tidak dapat diandalkan dan menunjukkan varians yang tinggi.

Salah satu cara untuk menghindari masalah ini adalah dengan menggunakan regresi komponen utama , yang menemukan M kombinasi linier (disebut “komponen utama”) dari prediktor p asli dan kemudian menggunakan kuadrat terkecil agar sesuai dengan model regresi linier yang menggunakan komponen utama sebagai prediktor.

Tutorial ini memberikan contoh langkah demi langkah tentang cara melakukan regresi komponen utama di R.

Langkah 1: Muat paket yang diperlukan

Cara termudah untuk melakukan regresi komponen utama di R adalah dengan menggunakan fungsi dalam paket pls .

 #install pls package (if not already installed)
install.packages(" pls ")

load pls package
library(pls)

Langkah 2: Sesuaikan model PCR

Untuk contoh ini, kita akan menggunakan kumpulan data R bawaan yang disebut mtcars yang berisi data tentang berbagai jenis mobil:

 #view first six rows of mtcars dataset
head(mtcars)

                   mpg cyl disp hp drat wt qsec vs am gear carb
Mazda RX4 21.0 6 160 110 3.90 2.620 16.46 0 1 4 4
Mazda RX4 Wag 21.0 6 160 110 3.90 2.875 17.02 0 1 4 4
Datsun 710 22.8 4 108 93 3.85 2.320 18.61 1 1 4 1
Hornet 4 Drive 21.4 6 258 110 3.08 3.215 19.44 1 0 3 1
Hornet Sportabout 18.7 8 360 175 3.15 3.440 17.02 0 0 3 2
Valiant 18.1 6 225 105 2.76 3,460 20.22 1 0 3 1

Untuk contoh ini, kami akan menyesuaikan model regresi komponen utama (PCR) dengan menggunakan hp sebagai variabel respons dan variabel berikut sebagai variabel prediktor:

• mpg
• menampilkan
• kotoran
• berat
• qdetik

Kode berikut menunjukkan cara menyesuaikan model PCR dengan data ini. Perhatikan argumen berikut:

• scale=TRUE : Ini memberi tahu R bahwa setiap variabel prediktor harus diskalakan agar memiliki rata-rata 0 dan deviasi standar 1. Hal ini memastikan bahwa tidak ada variabel prediktor yang memiliki pengaruh terlalu besar dalam model jika diukur dalam satuan yang berbeda. .
• validation=”CV” : Ini memberitahu R untuk menggunakan validasi k-fold cross untuk mengevaluasi kinerja model. Perhatikan bahwa ini menggunakan k=10 lipatan secara default. Perhatikan juga bahwa Anda dapat menentukan “LOOCV” untuk melakukan validasi silang Leave-One-Out .

 #make this example reproducible
set.seed(1)

#fit PCR model
model <- pcr(hp~mpg+disp+drat+wt+qsec, data=mtcars, scale= TRUE , validation=" CV ")

Langkah 3: Pilih jumlah komponen utama

Setelah kita menyesuaikan modelnya, kita perlu menentukan berapa banyak komponen utama yang layak dipertahankan.

Untuk melakukannya, cukup lihat uji root mean square error (uji RMSE) yang dihitung dengan validasi k-cross:

 #view summary of model fitting
summary(model)

Data: 
	Y dimension: 32 1
Fit method: svdpc
Number of components considered: 5

VALIDATION: RMSEP
Cross-validated using 10 random segments.
       (Intercept) 1 comp 2 comps 3 comps 4 comps 5 comps
CV 69.66 44.56 35.64 35.83 36.23 36.67
adjCV 69.66 44.44 35.27 35.43 35.80 36.20

TRAINING: % variance explained
    1 comp 2 comps 3 comps 4 comps 5 comps
X 69.83 89.35 95.88 98.96 100.00
hp 62.38 81.31 81.96 81.98 82.03

Hasilnya ada dua tabel menarik:

1. VALIDASI: RMSEP

Tabel ini menunjukkan pengujian RMSE yang dihitung dengan validasi k-fold cross. Kita dapat melihat hal berikut:

• Jika kita hanya menggunakan istilah asli dalam model, RMSE pengujiannya adalah 69.66 .
• Jika kita menambahkan komponen utama pertama, pengujian RMSE turun menjadi 44,56.
• Jika kita menambahkan komponen utama kedua, pengujian RMSE turun menjadi 35,64.

Kita dapat melihat bahwa menambahkan komponen utama tambahan sebenarnya menghasilkan peningkatan RMSE pengujian. Oleh karena itu, tampaknya optimal jika hanya menggunakan dua komponen utama dalam model akhir.

2. PELATIHAN: % varians dijelaskan

Tabel ini menunjukkan persentase variansi variabel respon yang dijelaskan oleh komponen utama. Kita dapat melihat hal berikut:

• Dengan hanya menggunakan komponen utama pertama, kita dapat menjelaskan 69,83% variasi variabel respon.
• Dengan menambahkan komponen utama kedua, kita dapat menjelaskan 89,35% variasi variabel respon.

Perhatikan bahwa kita masih dapat menjelaskan lebih banyak varians dengan menggunakan lebih banyak komponen utama, namun kita dapat melihat bahwa menambahkan lebih dari dua komponen utama tidak benar-benar meningkatkan persentase varians yang dijelaskan sebanyak itu.

Kita juga dapat memvisualisasikan pengujian RMSE (bersama dengan pengujian MSE dan R-squared) sebagai fungsi dari jumlah komponen utama menggunakan fungsi validationplot() .

 #visualize cross-validation plots
validationplot(model)
validationplot(model, val.type="MSEP")
validationplot(model, val.type="R2") 

Di setiap grafik, kita dapat melihat bahwa kecocokan model meningkat dengan menambahkan dua komponen utama, namun cenderung memburuk ketika kita menambahkan lebih banyak komponen utama.

Dengan demikian, model optimal hanya mencakup dua komponen utama pertama.

Langkah 4: Gunakan model akhir untuk membuat prediksi

Kita dapat menggunakan model PCR dua komponen utama terakhir untuk membuat prediksi tentang observasi baru.

Kode berikut menunjukkan cara membagi kumpulan data asli menjadi kumpulan pelatihan dan pengujian serta menggunakan model PCR dengan dua komponen utama untuk membuat prediksi pada kumpulan pengujian.

 #define training and testing sets
train <- mtcars[1:25, c("hp", "mpg", "disp", "drat", "wt", "qsec")]
y_test <- mtcars[26: nrow (mtcars), c("hp")]
test <- mtcars[26: nrow (mtcars), c("mpg", "disp", "drat", "wt", "qsec")]
    
#use model to make predictions on a test set
model <- pcr(hp~mpg+disp+drat+wt+qsec, data=train, scale= TRUE , validation=" CV ")
pcr_pred <- predict(model, test, ncomp= 2 )

#calculate RMSE
sqrt ( mean ((pcr_pred - y_test)^2))

[1] 56.86549

Kita melihat bahwa RMSE pengujian tersebut ternyata 56.86549 . Ini adalah deviasi rata-rata antara nilai hp yang diprediksi dan nilai hp yang diamati untuk observasi set pengujian.

Penggunaan penuh kode R dalam contoh ini dapat ditemukan di sini .