Bagaimana cara melakukan regresi ortogonal (kuadrat terkecil total) melalui PCA?

Saya selalu menggunakan lm()dalam R untuk melakukan regresi linier pada . Fungsi itu mengembalikan koefisien sehingga$y$$x$$\beta$

$$y = \beta x.$$

Hari ini saya belajar tentang kuadrat terkecil total dan princomp()fungsi itu (analisis komponen utama, PCA) dapat digunakan untuk melakukannya. Itu harus baik untuk saya (lebih akurat). Saya telah melakukan beberapa tes menggunakan princomp(), seperti:

r <- princomp( ~ x + y)

Masalah saya adalah: bagaimana menafsirkan hasilnya? Bagaimana saya bisa mendapatkan koefisien regresi? Dengan "koefisien" Maksud saya angka yang harus saya gunakan untuk mengalikan nilai untuk memberikan angka mendekati .x $\beta$$x$$y$

Kuadrat terkecil biasa vs kuadrat terkecil total

Pertama mari kita pertimbangkan kasus paling sederhana dari hanya satu variabel prediktor (independen) . Untuk kesederhanaan, biarkan keduanya x dan y dipusatkan, yaitu intersep selalu nol. Perbedaan antara regresi OLS standar dan regresi "orthogonal" TLS jelas ditunjukkan pada gambar (diadaptasi oleh saya) ini dari jawaban paling populer di utas terpopuler di PCA:$x$$x$$y$

OLS cocok persamaan dengan meminimalkan jarak kuadrat antara nilai-nilai yang diamati y dan nilai-nilai diprediksi y . TLS cocok dengan persamaan yang sama dengan meminimalkan jarak kuadrat antara ( x , y ) poin dan proyeksi mereka di telepon. Dalam kasus paling sederhana ini garis TLS hanyalah komponen utama pertama dari data 2D. Untuk menemukan β , jangan PCA pada ( x , y ) poin, yaitu membangun 2 × 2 kovarians matriks Σ dan menemukan vektor eigen pertama v$y=\beta x$$y$$\hat y$$(x,y)$$\beta$$(x,y)$$2\times 2$$\boldsymbol \Sigma$ ; maka β = v y / v x .$\mathbf v = (v_x, v_y)$$\beta = v_y/v_x$

Di Matlab:

 v = pca([x y]);    //# x and y are centered column vectors
 beta = v(2,1)/v(1,1);

Dalam R:

 v <- prcomp(cbind(x,y))$rotation
 beta <- v[2,1]/v[1,1]

By the way, ini akan menghasilkan kemiringan yang benar bahkan jika dan y tidak terpusat (karena fungsi PCA bawaan secara otomatis melakukan pemusatan). Untuk memulihkan intersep, hitung β 0 = ˉ y - β ˉ x .$x$$y$$\beta_0 = \bar y - \beta \bar x$

OLS vs. TLS, regresi berganda

Diberikan variabel dependen dan banyak variabel independen x i (sekali lagi, semuanya berpusat untuk kesederhanaan), regresi sesuai dengan persamaan y = β 1 x 1 + ... + β p x p . OLS tidak cocok dengan meminimalkan kesalahan kuadrat antara nilai-nilai yang diamati dari y dan diperkirakan nilai-nilai y . TLS melakukan fit dengan meminimalkan jarak kuadrat antara yang diamati ( x , y ) ∈ R p + $y$$x_i$

$$y= \beta_1 x_1 + \ldots + \beta_p x_p.$$ $y$$\hat y$$(\mathbf x, y)\in\mathbb R^{p+1}$ poin dan poin terdekat pada bidang regresi / pesawat terbang.

Perhatikan bahwa tidak ada lagi "garis regresi"! Persamaan di atas menentukan hyperplane : ini adalah bidang 2D jika ada dua prediktor, 3D hyperplane jika ada tiga prediktor, dll. Jadi solusi di atas tidak berfungsi: kita tidak bisa mendapatkan solusi TLS dengan hanya mengambil PC pertama (yang merupakan sebuah garis). Meski begitu, solusinya dapat dengan mudah diperoleh melalui PCA.

Seperti sebelumnya, PCA dilakukan pada titik . Ini menghasilkan p + 1 vektor eigen di kolom V . Vektor p eigen pertama mendefinisikan hyperplane p -dimensional H yang kita butuhkan; yang terakhir (nomor p + 1 ) vektor eigen v p + 1 adalah orthogonal untuk itu. Pertanyaannya adalah bagaimana mengubah basis H yang diberikan oleh vektor eigen p pertama menjadi koefisien β .$(\mathbf x, y)$$p+1$$\mathbf V$$p$$p$$\mathcal H$$p+1$$\mathbf v_{p+1}$$\mathcal H$$p$$\boldsymbol \beta$

Perhatikan bahwa jika kita mengatur untuk semua i ≠ k dan hanya x k = 1 , maka y = β k , yaitu vektor ( 0 , ... , 1 , ... , β k ) ∈ H terletak pada hyperplane H . Di sisi lain, kita tahu bahwa v p + 1 = ( v 1 , ... , v p + $x_i=0$$i \ne k$$x_k=1$$\hat y=\beta_k$

$$(0,\ldots, 1, \ldots, \beta_k) \in \mathcal H$$ $\mathcal H$ ortogonal untuk itu. Yaitu produk titik mereka harus nol: v k + β k v p + 1 = 0 ⇒ β k = - v k / v p + 1 $$\mathbf v_{p+1}=(v_1, \ldots, v_{p+1}) \:\bot\: \mathcal H$$ $$v_k + \beta_k v_{p+1}=0 \Rightarrow \beta_k = -v_k/v_{p+1}.$$

Di Matlab:

 v = pca([X y]);    //# X is a centered n-times-p matrix, y is n-times-1 column vector
 beta = -v(1:end-1,end)/v(end,end);

Dalam R:

 v <- prcomp(cbind(X,y))$rotation
 beta <- -v[-ncol(v),ncol(v)] / v[ncol(v),ncol(v)]

Sekali lagi, ini akan menghasilkan lereng yang benar bahkan jika dan y tidak terpusat (karena fungsi PCA bawaan secara otomatis melakukan pemusatan). Untuk memulihkan intersep, hitung β 0 = ˉ y - ˉ x β .$x$$y$$\beta_0 = \bar y - \bar {\mathbf x} \boldsymbol \beta$

$x$$(x,y)$$v^{(1)}_y/v^{(1)}_x=-v^{(2)}_x/v^{(2)}_y$

Solusi formulir tertutup untuk TLS

$\boldsymbol \beta$

$\mathbf X$$\mathbf y$$\mathbf v_{p+1}$$[\mathbf X\: \mathbf y]$$\sigma^2_{p+1}$$-\mathbf v_{p+1}/v_{p+1} = (\boldsymbol \beta\:\: -1)^\top$

$$\left(\begin{array}{c}\mathbf X^\top \mathbf X & \mathbf X^\top \mathbf y\\ \mathbf y^\top \mathbf X & \mathbf y^\top \mathbf y\end{array}\right) \left(\begin{array}{c}\boldsymbol \beta \\ -1\end{array}\right) = \sigma^2_{p+1}\left(\begin{array}{c}\boldsymbol \beta \\ -1\end{array}\right),$$ $$\boldsymbol \beta_\mathrm{TLS} = (\mathbf X^\top \mathbf X - \sigma^2_{p+1}\mathbf I)^{-1} \mathbf X^\top \mathbf y,$$ $$\boldsymbol \beta_\mathrm{OLS} = (\mathbf X^\top \mathbf X)^{-1} \mathbf X^\top \mathbf y.$$

Regresi berganda multivarian

Rumus yang sama dapat digeneralisasi ke kasus multivariat, tetapi bahkan untuk menentukan apa yang dilakukan TLS multivariat, akan memerlukan beberapa aljabar. Lihat Wikipedia di TLS . Regresi OLS multivariat setara dengan sekelompok regresi OLS univariat untuk setiap variabel dependen, tetapi dalam kasus TLS tidak demikian.

Berdasarkan implementasi GNU Octave yang naif ditemukan di sini , sesuatu seperti ini mungkin (butiran garam, sudah terlambat) bekerja.

tls <- function(A, b){

  n <- ncol(A)
  C <- cbind(A, b)

  V <- svd(C)$v
  VAB <- V[1:n, (n+1):ncol(V)]
  VBB <- V[(n+1):nrow(V), (n+1):ncol(V)]
  return(-VAB/VBB)
}

princompmenjalankan analisis komponen utama alih-alih total kuadrat terkecil regresi. Sejauh yang saya tahu tidak ada fungsi R atau paket yang melakukan TLS; paling banyak ada regresi Deming di MethComp .
Namun, tolong perlakukan ini sebagai saran bahwa kemungkinan besar itu tidak sepadan.