Komputasi Regresi

Alfa Nugraha¹

2023-11-14

Regresi Linear dalam bentuk matriks

Model regresi linear dapat dituliskan sebagai berikut:

\[ y_{i} = \beta_{0}+\beta_{1}x_{1}+\beta_{2}x_{2}+\ldots + \beta_{p-1}x_{p-1}+\epsilon_{i} \] dengan

\(y_{i}=\) peubah respon

\(x_{1},x_{2},\ldots,x_{p-1}=\) peubah explanatory (penjelas)

\(\beta_{0},\beta_{1},\ldots,\beta_{p-1}=\) koefisien peubah explanatory (penjelas)

\(\epsilon_{i}=\) galat

Model regresi linear ini dapat dituliskan dalam bentuk matriks sebagai berikut

\[ y_{n\times1}=X_{n\times p}\beta_{p\times 1}+ \epsilon_{n\times1} \] dengan \(X_{n\times p}\) merupakan design matrix atau model matrix

\[ X_{n\times p} = \begin{pmatrix} 1 & X_{1,1} & X_{1,2} & \ldots & X_{1,p-1} \\ 1 & X_{2,1} & X_{2,2}& \ldots & X_{2,p-1} \\ \vdots & \vdots & \vdots &\ddots & \vdots \\ 1 & X_{n,1} & X_{n,2}& \ldots & X_{n,p-1} \end{pmatrix} \]

\(\beta_{p\times 1}\) merupakan vektor peubah respon

\[ \beta_{p\times 1} = \begin{pmatrix} \beta_{0} \\ \beta_{1} \\ \beta_{2} \\ \vdots \\ \beta_{p-1} \end{pmatrix} \]

\(y_{n\times 1}\) merupakan vektor peubah respon

\[ y_{n\times 1} = \begin{pmatrix} y_{1} \\ y_{2} \\ \vdots \\ y_{n} \end{pmatrix} \]

\(\epsilon_{n\times1}\) merupakan vektor peubah galat

\[ \epsilon_{n\times1} = \begin{pmatrix} \epsilon_{1} \\ \epsilon_{2} \\ \vdots \\ \epsilon_{n} \end{pmatrix} \]

Pendugaan parameter koefisien regresi dapat dilakukan dengan menggunakan metode kuadrat terkecil atau least square. Hasil pendugaan koefisien \(b\) didapatkan dengan meminimumkan Jumlah Kuardat Galat (JKG) atau Sum of Squared Error (SSE)

\[ \text{SSE}=(y-Xb)'(y-Xb) \] sehingga diperoleh persamaan normal (normal equation)

\[ (X'X)b= X'y \]

Solusi dari persamaan normal tersebut yang akan menjadi dugaan koefisien regresi \(b\)

\[ b = (X'X)^{-1}X'y \]

Matriks Simetrik dan Definit Positive

Matriks Simetrik

\(A\) merupakan matriks berukuran \(n\times n\). \(A\) dikatakan simetrik jika dan hanya jika

\[ A=A' \space \text{dengan} \space a_{ij}=a_{ji} \]

Matriks Definit positif

\(A\) merupakan matriks berukuran \(n\times n\) dan

\[ {x} = \begin{pmatrix} x_{1} \\ x_{2} \\ \vdots \\ x_{n} \end{pmatrix} \] merupakan vektor dengan ukuruan \(n\times 1\). Maka \(q=x'Ax\) disebut bentuk kuadrat dalam \(x\). Kemudian, \(A\) disebut matriks bentuk kuadrat (quadratic form)

Bentuk kuadrat \(x'Ax\) dikatakan definit positif jika

\[ x'Ax>0, \space \text{untuk setiap} \space x\neq 0 \]

dan \(x'Ax\) dikatakan semidefinit positif jika

\[ x'Ax\geq0, \space \text{untuk setiap} \space x\neq 0 \] dan

\[ x'Ax=0, \space \text{untuk suatu} \space x\neq 0 \]

Istilah ini juga berlaku untuk matriks. Jika suatu matriks definit positif maka bentuk kuadrat-nya juga definit positif. Begitupun jika Jika suatu matriks semidefinit positif maka bentuk kuadrat-nya juga semidefinit positif

SWEEP Operator

SWEEP operator melakukan operasi baris elementer pada sistem persamaan linier. SWEEP operator digunakan untuk membuat model regresi dengan melakukan sweep baris tertentu dari design matrix \(X'X\). SWEEP operator juga dipakai untuk mendapatkan dugaan koefisien regresi, jumlah kuadrat galat (JKG) atau Sum of Squares Error (SSE), dan generalized inverse secara bersamaan. Hal ini berlaku tidak hanya untuk variabel respon tunggal tetapi juga untuk beberapa respon.

SWEEP operator mirip seperti eliminasi Gauss yang merupakan teknik umum untuk menyelesaikan sistem persamaan linear. SWEEP operator paling sering digunakan untuk menyelesaikan masalah kuadrat terkecil atau least square untuk persamaan \(X'X = X'y\). Oleh karena itu, dalam praktiknya, SWEEP operator biasanya diterapkan pada baris-baris matriks definit positif simetris.

Algoritma SWEEP Operator

Goodnight (1979) mendefinisikan SWEEP Operator sebagai urutan operasi baris. Diberikan matriks definit positif simetris \(A\), maka \(\text{SWEEP}(A, k)\) memodifikasi matriks \(A\) dengan menggunakan elemen pivot \(A[k,k]\) dan baris ke-\(k\), sebagai berikut:

Tetapkan \(D = A[k,k]\), dimana \(A[k,k]\) adalah elemen diagonal ke-\(k\).
Bagilah baris ke-\(k\) dengan \(D\).
Untuk setiap \(i \neq k\), lakukan:
1. Tetapkan \(B = A[i,k]\)
2. Kurangkan baris \(i\) dengan \(B \times A[k,]\)
3. Tetapkan \(A[i,k] = –B/D\).
Tetapkan \(A[k,k] = 1/D\).

Matriks \(A\) didefinisikan dalam matriks partisi (partitioned Matrix/block Matrix)

\[ A = \begin{pmatrix} X'X & |& X'y \\ y'X & |& y'y \end{pmatrix} \]

dimana \(A\) bisa kita hitung dengan

\[ A = M'M \]

Matriks \(M\) merupakan matrik gabungan antara design matrix atau model matrix \(X\) dengan vektor respon \(y\)

\[ M = (X|y) \]

dimana matrix \(X\) didefinisikan

dan

\[ y_{n\times 1} = \begin{pmatrix} y_{1} \\ y_{2} \\ \vdots \\ y_{n} \end{pmatrix} \]

sehingga matrix \(M\) bisa ditulis

\[ M_{n\times (p+1)} = \begin{pmatrix} 1 & X_{1,1} & X_{1,2} & \ldots & X_{1,p-1} & y_{1}\\ 1 & X_{2,1} & X_{2,2}& \ldots & X_{2,p-1} & y_{2}\\ \vdots & \vdots & \vdots &\ddots & \vdots & \vdots\\ 1 & X_{n,1} & X_{n,2}& \ldots & X_{n,p-1} & y_{n} \end{pmatrix} \]

Hasil akhir dari SWEEP OPERATOR adalah matriks yang didefinisikan

\[ \text{SWEEP}(A,k) = \begin{pmatrix} (X'X)^{-1} & |& (X'X)^{-1}X'y \\ -y'X (X'X)^{-1}& |& y'y-y'X (X'X)^{-1}X'y \end{pmatrix} \]

atau bisa ditulis

\[ \text{SWEEP}(A,k) = \begin{pmatrix} (X'X)_{p\times p}^{-1} & |& b_{p\times 1} \\ -b_{1\times p}'& |& \text{SSE}_{1\times1} \end{pmatrix} \]

Setiap langkah pada SWEEP Operator akan mendapatkan nilai koefisien untuk peubah penjelas dan nilai JKG. Hal ini bisa dirangkum pada tabel berikut:

Tahap	Persamaan	Koefisien	JKG
1	\(y=b_{0}\)	\(b_{0}\)	\(\text{JKG}_{0}\)
2	\(y=b_{0}+b_{1}x_{1}\)	\(b_{0},b_{1}\)	\(\text{JKG}_{1}\)
3	\(y=b_{0}+b_{1}x_{1}+b_{2}x_{2}\)	\(b_{0},b_{1},b_{2}\)	\(\text{JKG}_{2}\)
\(\vdots\)	\(\vdots\)	\(\vdots\)	\(\vdots\)
k	\(y=b_{0}+b_{1}x_{1}+b_{2}x_{2}+\ldots+b_{p}x_{p}\)	\(b_{0},b_{1},b_{2},\ldots,b_{p}\)	\(\text{JKG}_{k-1}\)

Untuk menjalankan SWEEP Operator diatas dapat menggunakan fungsi sweep.operator dari package fastmatrix

Dekomposisi Cholesky

Jika \(A\) merupakan matriks simetrik definit positif berukuran \(n\times n\), maka terdapat suatu matriks segitiga atas \(L\) dengan diagonal utama positif sedemikian sehingga

\[ A=LL' \]

dengan

\[ L=\begin{pmatrix} l_{1,1} & 0 & \ldots & 0 \\ l_{2,1} & l_{2,2}& \ldots & 0 \\ \vdots & \vdots &\ddots & \vdots \\ l_{n,1} & l_{n,2}& \ldots & l_{n,n} \end{pmatrix} \]

Dekomposisi matriks \(A\) menjadi \(LL'\) dinamakan Dekomposisi Cholesky. Matriks \(L\) dapat dihitung dengan formula berikut:

Untuk baris ke\(k\), dapat dihitung

\[ l_{ki}=\frac{a_{ki}-\Sigma_{j=1}^{i-1}l_{ij}l_{kj}}{l_{ii}}, \space \text{untuk} \space i=1,2,\ldots,k-1 \]

dan

\[ l_{kk} = \sqrt{a_{kk}-\Sigma_{j=1}^{k-1}l_{kj}^2} \]

Dekomposisi Cholesky dapat digunakan untuk menghitung koefisien regresi dengan menyelesaikan persamaan normal (normal equation) berikut

\[ (X'X)b= X'y \]

karena matriks \(X'X\) merupakan matriks simetrik definit positif berukuran \(p\times p\) maka

\[ X'X=LL' \]

sehingga

\[ (LL')b= X'y \]

Koefisien regresi \(b\) bisa diperoleh secara langsung

\[ b=L'^{-1}L^{-1}X'y \]

atau dengan memisalkan \(L'b=z\) dan menerapkan dua langkah berikut:

Hitung solusi untuk \(z\) dengan

\[ z=L^{-1}X'y \] Langkah 1 ini dinamakan forward solve lower triangular system

Berdasarkankan hasil yang diperoleh pada 1, hitung

\[ b=L'^{-1}z \] Langkah 2 ini dinamakan backward solve upper triangular system

Apakah ada manfaat menggunakan matriks segitiga atas atau bawah dalam menyelesaikan sistem persamaan linear (persamaan normal regresi)? Ya. Terdapat fungsi di R yang cepat waktu komputasi-nya untuk menyelesaikan sistem persamaan linier ketika matriks dalam persamaan adalah marix segitiga atas atau bawah, yaitu fungsi backsolve dan forwardsolve. Berikut adalah user-defined function menggunakan fungsi backsolve dan forwardsolve

lm_chol <- function(X, y){
  
    # menghitung X'X
    XX <- crossprod(X)
    # menghitung X'Y
    Xy <- crossprod(X, y)
    
    # cholesky decomposition
    upper <- chol(XX)
    
    #menghitung forward lower triangular system
    #dan backward solve upper triangular system
    
    # Kita perlu melakukan hal berikut
    # z <- forwardsolve( t(upper), XY )
    # namun  transpose dari L bisa dihindari dengan
    # argumnen transpose=TRUE pada fungsi backsolve
    z <- backsolve(upper, Xy, transpose = TRUE)
    coef <- as.vector(backsolve(upper, z))
    
    return(coef)
}

Soal dan Pembahasan

Soal 1

Berdasarkan suatu data yang ingin dimodelkan dengan regresi linear, diketahui informasi-informasi berikut

\[ X'X= \begin{pmatrix} 6 & 0 & 12 \\ 0 & 6 & 0 \\ 12 & 0 & 28 \end{pmatrix} \]

\[ X'y= \begin{pmatrix} 12 \\ 2 \\ 25 \end{pmatrix} \]

\[ y'y= 28 \]

Hitunglah penduga koefisien regresi dengan menggunakan SWEEP Operator tanpa menggunakan bantuan software!
Hitunglah penduga koefisien regresi dengan menggunakan SWEEP Operator dengan menggunakan sweep.operator di R!
Buatlah tabel Tahapan SWEEP operator!

Jawaban

Hitunglah penduga koefisien regresi dengan menggunakan SWEEP Operator tanpa menggunakan bantuan software!

Langkah Pertama: Membuat Matriks A

\[ A = \begin{pmatrix} X'X & |& X'y \\ y'X & |& y'y \end{pmatrix} \]

\[ A = \begin{pmatrix} 6 & 0 & 12 & 12\\ 0 & 6 & 0 & 2\\ 12 & 0 & 28 & 25 \\ 12 & 2 & 25 & 28 \end{pmatrix} \]

Langkah Kedua: Lakukan Tahapan-tahapan pada SWEEP operator

Tahapan Berikut untuk \(\text{SWEEP}(A,k=1)\). Misalkan saja \(A^{*}=\text{SWEEP}(A,k=1)\)

Tahap 1

Menentukan \(D=A^{*}[1,1]\), sehingga \(D=6\)

Tahap 2

Membagi baris ke-\(k=1\) dengan \(D=6\), sehingga

\[ A^{*} = \begin{pmatrix} \boldsymbol{6/6} & \boldsymbol{0/6} & \boldsymbol{12/6} & \boldsymbol{12/6}\\ 0 & 6 & 0 & 2\\ 12 & 0 & 28 & 25 \\ 12 & 2 & 25 & 28 \end{pmatrix} \]

\[ A^{*} = \begin{pmatrix} \boldsymbol{1} & \boldsymbol{0} & \boldsymbol{2} & \boldsymbol{2}\\ 0 & 6 & 0 & 2\\ 12 & 0 & 28 & 25 \\ 12 & 2 & 25 & 28 \end{pmatrix} \]

Tahap 3

Untuk Baris \(A^{*}[i=2,]\)

Tetapkan \(B=A^{*}[2,1]=0\)
Kurangkan baris ke-\(i=2\) dengan \(B \times A^{*}[k=1,]\)

\[A^{*}[2,]-B \times A^{*}[1,]\]

\[\begin{pmatrix} 0 & 6 & 0 & 2\\ \end{pmatrix}-0\times \begin{pmatrix} \boldsymbol{1} & \boldsymbol{0} & \boldsymbol{2} & \boldsymbol{2}\\ \end{pmatrix}\]

\[\begin{pmatrix} 0 & 6 & 0 & 2\\ \end{pmatrix}\]

sehingga

\[ A^{*} = \begin{pmatrix} 1 & 0 & 2 & 2\\ \boldsymbol{0} & \boldsymbol{6} & \boldsymbol{0} & \boldsymbol{2}\\ 12 & 0 & 28 & 25 \\ 12 & 2 & 25 & 28 \end{pmatrix} \]

Tetapkan \(A^{*}[i=2,k=1] = –B/D\).

\[ A^{*}[2,1]=-0/6=0 \]

sehingga

\[ A^{*} = \begin{pmatrix} 1 & 0 & 2 & 2\\ \boldsymbol{0} & \boldsymbol{6} & \boldsymbol{0} & \boldsymbol{2}\\\ 12 & 0 & 28 & 25 \\ 12 & 2 & 25 & 28 \end{pmatrix} \]

Untuk Baris \(A^{*}[i=3,]\)

Tetapkan \(B=A^{*}[3,1]=12\)
Kurangkan baris ke-\(i=3\) dengan \(B \times A^{*}[k=1,]\)

\[A^{*}[3,]-B \times A^{*}[1,]\]

\[\begin{pmatrix} 12 & 0 & 28 & 25\\ \end{pmatrix}-12\times \begin{pmatrix} \boldsymbol{1} & \boldsymbol{0} & \boldsymbol{2} & \boldsymbol{2}\\ \end{pmatrix}\]

\[\begin{pmatrix} 0 & 0 & 4 & 1\\ \end{pmatrix}\]

sehingga

\[ A^{*} = \begin{pmatrix} 1 & 0 & 2 & 2\\ 0 & 6 & 0 & 2\\ \boldsymbol{0} & \boldsymbol{0} & \boldsymbol{4} & \boldsymbol{1} \\ 12 & 2 & 25 & 28 \end{pmatrix} \]

Tetapkan \(A^{*}[i=3,k=1] = –B/D\).

\[ A^{*}[3,1]=-12/6=-2 \]

sehingga

\[ A^{*} = \begin{pmatrix} 1 & 0 & 2 & 2\\ 0 & 6 & 0 & 2\\ \boldsymbol{-2} & \boldsymbol{0} & \boldsymbol{4} & \boldsymbol{1} \\ 12 & 2 & 25 & 28 \end{pmatrix} \]

Untuk Baris \(A^{*}[i=4,]\)

Tetapkan \(B=A^{*}[3,1]=12\)
Kurangkan baris ke-\(i=4\) dengan \(B \times A^{*}[k=1,]\)

\[A^{*}[4,]-B \times A^{*}[1,]\]

\[\begin{pmatrix} 12 & 2 & 25 & 28\\ \end{pmatrix}-12\times \begin{pmatrix} \boldsymbol{1} & \boldsymbol{0} & \boldsymbol{2} & \boldsymbol{2}\\ \end{pmatrix}\]

\[\begin{pmatrix} 0 & 2 & 1 & 4\\ \end{pmatrix}\]

sehingga

\[ A^{*} = \begin{pmatrix} 1 & 0 & 2 & 2\\ 0 & 6 & 0 & 2\\ -2 & 0 & 4 & 1 \\ \boldsymbol{0} & \boldsymbol{2} & \boldsymbol{1} & \boldsymbol{4} \end{pmatrix} \]

Tetapkan \(A^{*}[i=4,k=1] = –B/D\).

\[ A^{*}[4,1]=-12/6=-2 \]

sehingga

\[ A^{*} = \begin{pmatrix} 1 & 0 & 2 & 2\\ 0 & 6 & 0 & 2\\ -2 & 0 & 4 & 1 \\ \boldsymbol{-2} & \boldsymbol{2} & \boldsymbol{1} & \boldsymbol{4} \end{pmatrix} \]

Tahap 4

Tetapkan \(A^{*}[k,k] = 1/D=1/6=0.1667\)

Sehingga diperoleh

\[ A^{*} = \begin{pmatrix} \boldsymbol{0.1667} & 0 & 2 & 2\\ 0 & 6 & 0 & 2\\ -2 & 0 & 4 & 1 \\ -2 & 2 & 1 & 4 \end{pmatrix} \]

Jadi Matriks Akhir yang diperoleh

\[ \text{SWEEP}(A,k=1)=A^{*} = \begin{pmatrix} 0.1667 & 0 & 2 & 2\\ 0 & 6 & 0 & 2\\ -2 & 0 & 4 & 1 \\ -2 & 2 & 1 & 4 \end{pmatrix} \]

intercept regresi bisa diambil dari kolom ke 4 baris 1. Hal ini sesuai dengan definisi matriks hasil akhir dari SWEEP Operator

\[ \text{SWEEP}(A,k) = \begin{pmatrix} (X'X)_{p\times p}^{-1} & |& b_{p\times 1} \\ -b_{1\times p}'& |& \text{SSE}_{1\times1} \end{pmatrix} \]

Sehingga persamaan regresi yang diperoleh \(y=2\) dengan \(SSE=4\)

Tahapan Berikut untuk \(\text{SWEEP}(A,k=2)\). Misalkan saja \(A^{*}=\text{SWEEP}(A,k=1)\) dan

\(A^{**}=\text{SWEEP}(A^{*},k=2)\)

Tahap 1

Menentukan \(D=A^{**}[2,2]\), sehingga \(D=6\)

Tahap 2

Membagi baris ke-\(k=2\) dengan \(D=6\), sehingga

\[ A^{**} = \begin{pmatrix} 0.1667 & 0 & 2 & 2\\ \boldsymbol{0/6} & \boldsymbol{6/6} & \boldsymbol{0/6} & \boldsymbol{2/6}\\ -2 & 0 & 4 & 1 \\ -2 & 2 & 1 & 4 \end{pmatrix} \]

\[ A^{**} = \begin{pmatrix} 0.1667 & 0 & 2 & 2\\ \boldsymbol{0} & \boldsymbol{1} & \boldsymbol{0} & \boldsymbol{0.333}\\ -2 & 0 & 4 & 1 \\ -2 & 2 & 1 & 4 \end{pmatrix} \]

Tahap 3

Untuk Baris \(A^{**}[i=1,]\)

Tetapkan \(B=A^{**}[1,2]=0\)
Kurangkan baris ke-\(i=1\) dengan \(B \times A^{**}[k=2,]\)

\[A^{**}[1,]-B \times A^{*}[2,]\]

\[\begin{pmatrix} 0.1667 & 0 & 2 & 2\\ \end{pmatrix}-0\times \begin{pmatrix} \boldsymbol{0} & \boldsymbol{1} & \boldsymbol{0} & \boldsymbol{0.333}\\ \end{pmatrix}\]

\[\begin{pmatrix} 0.1667 & 0 & 2 & 2\\ \end{pmatrix}\]

sehingga

\[ A^{**} = \begin{pmatrix} \boldsymbol{0.1667} & \boldsymbol{0} & \boldsymbol{2} & \boldsymbol{2}\\ 0 & 1 & 0 & 0.333\\ -2 & 0 & 4 & 1 \\ -2 & 2 & 1 & 4 \end{pmatrix} \]

Tetapkan \(A^{**}[i=1,k=2] = –B/D\).

\[ A^{**}[1,2]=-0/6=0 \]

sehingga

\[ A^{**} = \begin{pmatrix} \boldsymbol{0.1667} & \boldsymbol{0} & \boldsymbol{2} & \boldsymbol{2}\\ 0 & 1 & 0 & 0.333\\ -2 & 0 & 4 & 1 \\ -2 & 2 & 1 & 4 \end{pmatrix} \]

Untuk Baris \(A^{**}[i=3,]\)

Tetapkan \(B=A^{**}[3,2]=0\)
Kurangkan baris ke-\(i=3\) dengan \(B \times A^{**}[k=2,]\)

\[A^{**}[3,]-B \times A^{**}[2,]\]

\[\begin{pmatrix} -2 & 0 & 4 & 1 \\ \end{pmatrix}-0\times \begin{pmatrix} \boldsymbol{0} & \boldsymbol{1} & \boldsymbol{0} & \boldsymbol{0.333}\\ \end{pmatrix}\]

\[\begin{pmatrix} -2 & 0 & 4 & 1 \\ \end{pmatrix}\]

sehingga

\[ A^{*} = \begin{pmatrix} 0.1667 & 0 & 2 & 2\\ 0 & 1 & 0 & 0.333\\ \boldsymbol{-2} & \boldsymbol{0} & \boldsymbol{4} & \boldsymbol{1} \\ -2 & 2 & 1 & 4 \end{pmatrix} \]

Tetapkan \(A^{**}[i=3,k=2] = –B/D\).

\[ A^{**}[3,2]=-0/6=0 \]

sehingga

\[ A^{**} = \begin{pmatrix} 0.1667 & 0 & 2 & 2\\ 0 & 1 & 0 & 0.333\\ \boldsymbol{-2} & \boldsymbol{0} & \boldsymbol{4} & \boldsymbol{1} \\ -2 & 2 & 1 & 4 \end{pmatrix} \]

Untuk Baris \(A^{**}[i=4,]\)

Tetapkan \(B=A^{*}[4,1]=2\)
Kurangkan baris ke-\(i=4\) dengan \(B \times A^{**}[k=2,]\)

\[A^{**}[4,]-B \times A^{**}[2,]\]

\[\begin{pmatrix} -2 & 2 & 1 & 4\\ \end{pmatrix}-2\times \begin{pmatrix} \boldsymbol{0} & \boldsymbol{1} & \boldsymbol{0} & \boldsymbol{0.333}\\ \end{pmatrix}\]

\[\begin{pmatrix} -2 & 0 & 1 & 3.333\\ \end{pmatrix}\]

sehingga

\[ A^{**} = \begin{pmatrix} 0.1667 & 0 & 2 & 2\\ 0 & 1 & 0 & 0.333\\ -2 & 0 & 4 & 1 \\ \boldsymbol{-2} & \boldsymbol{0} & \boldsymbol{1} & \boldsymbol{3.33} \end{pmatrix} \]

Tetapkan \(A^{**}[i=4,k=2] = –B/D\).

\[ A^{**}[4,1]=-2/6=-0.333 \]

sehingga

\[ A^{**} = \begin{pmatrix} 0.1667 & 0 & 2 & 2\\ 0 & 1 & 0 & 0.333\\ -2 & 0 & 4 & 1 \\ \boldsymbol{-2} & \boldsymbol{-0.333} & \boldsymbol{1} & \boldsymbol{3.33} \end{pmatrix} \]

Tahap 4

Tetapkan \(A^{**}[k,k] = 1/D=1/6=1.667\)

Sehingga diperoleh

\[ A^{**} = \begin{pmatrix} 0.1667 & 0 & 2 & 2\\ 0 & \boldsymbol{0.1667}& 0 & 0.333\\ -2 & 0 & 4 & 1 \\ -2 & -0.333 & 1 & 3.333 \end{pmatrix} \]

Jadi Matriks Akhir yang diperoleh

\[ \text{SWEEP}(A^{*},k=2)=A^{**} = \begin{pmatrix} 0.1667 & 0 & 2 & 2\\ 0 & 0.1667 & 0 & 0.333\\ -2 & 0 & 4 & 1 \\ -2 & -0.333 & 1 & 3.333 \end{pmatrix} \]

intercept dan koefisien regresi bisa diambil dari kolom ke 4 baris 2. Hal ini sesuai dengan definisi matriks hasil akhir dari SWEEP Operator

\[ \text{SWEEP}(A,k) = \begin{pmatrix} (X'X)_{p\times p}^{-1} & |& b_{p\times 1} \\ -b_{1\times p}'& |& \text{SSE}_{1\times1} \end{pmatrix} \]

Sehingga persamaan regresi yang diperoleh \(y=2+0.333x_{1}\) dengan \(SSE=3.333\)

Tahapan Berikut untuk \(\text{SWEEP}(A,k=3)\). Silahkan lanjutkan sebagai latihan

Hitunglah penduga koefisien regresi dengan menggunakan SWEEP Operator dengan menggunakan sweep.operator di R!

Langkah pertama : Input Matriks X’X, X’y dan y’y ke R

XX = matrix(c(6,0,12,
              0,6,0,
              12,0,28),
            nrow = 3,byrow = TRUE)
XX

##      [,1] [,2] [,3]
## [1,]    6    0   12
## [2,]    0    6    0
## [3,]   12    0   28

Xy <- matrix(c(12,
               2,
               25),
             ncol = 1)
Xy

##      [,1]
## [1,]   12
## [2,]    2
## [3,]   25

yy <- 28

Langkah kedua: Membuat Matriks A

A <- rbind(cbind(XX,Xy),cbind(t(Xy),yy))
colnames(A)

## [1] ""   ""   ""   "yy"

# hapus nama kolom supaya menjadi 
# matriks simetris
colnames(A) <- NULL
A

##      [,1] [,2] [,3] [,4]
## [1,]    6    0   12   12
## [2,]    0    6    0    2
## [3,]   12    0   28   25
## [4,]   12    2   25   28

Langkah Ketiga : Hitung koefisien Regresi dengan fungsi sweep.operator

Karena ukuran matriks \(X'X\) adalah \(3 \times 3\) maka banyaknya koefisien regresi yang bisa diduga ada 3 dengan 2 koefisien untuk peubah penjelas dan satu untuk intercept

res_sweep <- sweep.operator(A,k = 1:3)
dim(res_sweep)

## [1] 4 4

koefisien regresi bisa diambil dari kolom ke 4 baris 1 sampai 4. Hal ini sesuai dengan definisi matriks hasil akhir dari SWEEP Operator

\[ \text{SWEEP}(A,k) = \begin{pmatrix} (X'X)_{p\times p}^{-1} & |& b_{p\times 1} \\ -b_{1\times p}'& |& \text{SSE}_{1\times1} \end{pmatrix} \]

coef_sweep <- res_sweep[-4,4]
length(coef_sweep)

## [1] 3

coef_sweep

## [1] 1.5000000 0.3333333 0.2500000

Buatlah tabel Tahapan SWEEP operator!

Untuk membuat tabel Tahapan SWEEP Operator, kita perlu mengatur nilai k

Jika k=1 maka kita hanya mendapatkan dugaan intercept b0

coef_sweep_b0 <- sweep.operator(A,k = 1)
coef_sweep_b0

##            [,1] [,2] [,3] [,4]
## [1,]  0.1666667    0    2    2
## [2,]  0.0000000    6    0    2
## [3,] -2.0000000    0    4    1
## [4,] -2.0000000    2    1    4

## intercept
coef_sweep_b0[1,4]

## [1] 2

## SEE
coef_sweep_b0[4,4]

## [1] 4

Jika k=1:2 maka kita hanya mendapatkan dugaan intercept b0 dan koefisien b1

coef_sweep_b1 <- sweep.operator(A,k = 1:2)
coef_sweep_b1

##            [,1]       [,2] [,3]      [,4]
## [1,]  0.1666667  0.0000000    2 2.0000000
## [2,]  0.0000000  0.1666667    0 0.3333333
## [3,] -2.0000000  0.0000000    4 1.0000000
## [4,] -2.0000000 -0.3333333    1 3.3333333

## intercept+b1
coef_sweep_b1[1:2,4]

## [1] 2.0000000 0.3333333

## SEE
coef_sweep_b1[4,4]

## [1] 3.333333

Jika k=1:3 maka kita hanya mendapatkan dugaan intercept b0, koefisien b1 dan b2

coef_sweep_b2 <- sweep.operator(A,k = 1:3)
coef_sweep_b2

##           [,1]       [,2]  [,3]      [,4]
## [1,]  1.166667  0.0000000 -0.50 1.5000000
## [2,]  0.000000  0.1666667  0.00 0.3333333
## [3,] -0.500000  0.0000000  0.25 0.2500000
## [4,] -1.500000 -0.3333333 -0.25 3.0833333

## intercept+b1+b2
coef_sweep_b2[1:3,4]

## [1] 1.5000000 0.3333333 0.2500000

## SEE
coef_sweep_b2[4,4]

## [1] 3.083333

Sehingga kita bisa mengisi tabel tahapan SWEEP operator

Tahap	Persamaan	Koefisien	JKG
1	\(y=2\)	\(2\)	\(4\)
2	\(y=2+0.33x_{1}\)	\(2,0.33\)	\(3.33\)
3	\(y=1.5+0.33x_{1}+0.25x_{2}\)	\(1.5,0.33,0.25\)	\(3.083\)

Soal 2

Berdasarkan suatu data yang ingin dimodelkan dengan regresi linear, diketahui informasi-informasi berikut

\[ X'X= \begin{pmatrix} 32 & 198 & 115.09 \\ 198 & 1324 & 691.4 \\ 115.09 & 691.4 & 422.7907 \end{pmatrix} \]

\[ X'y= \begin{pmatrix} 642.9 \\ 3693.6 \\ 2380.277 \end{pmatrix} \]

\[ y'y= 14042 \]

Hitunglah penduga koefisien regresi dengan menggunakan Dekomposisi Cholesky dengan menggunakan R!

Jawaban

Hitunglah penduga koefisien regresi dengan menggunakan Dekomposisi Cholesky dengan menggunakan R!

Langkah pertama : Input Matriks X’X dan X’y ke R

XX = matrix(c(32,198,115.09,
              198,1324,691.4,
              115.09, 691.4, 422.7907), nrow = 3,byrow = TRUE)
XX

##        [,1]   [,2]     [,3]
## [1,]  32.00  198.0 115.0900
## [2,] 198.00 1324.0 691.4000
## [3,] 115.09  691.4 422.7907

Xy <- matrix(c(642.9,
               3693.6,
               2380.277), ncol = 1)
Xy

##          [,1]
## [1,]  642.900
## [2,] 3693.600
## [3,] 2380.277

Langkah kedua: menghitung dekomposisi Cholesky

upper <- chol(XX)
upper

##          [,1]      [,2]      [,3]
## [1,] 5.656854 35.001786 20.345230
## [2,] 0.000000  9.943591 -2.083691
## [3,] 0.000000  0.000000  2.126159

forward lower triangular system dan backward solve upper triangular system

# Kita perlu melakukan hal berikut
# z <- forwardsolve( t(upper), XY )
# namun  transpose dari L bisa dihindari dengan
# argumnen transpose=TRUE pada fungsi backsolve
z <- backsolve(upper, Xy, transpose = TRUE)
coef <- as.vector(backsolve(upper, z))
coef

## [1] 28.724665 -2.483514  1.871983

Soal 3

Dengan menggunakan data longley dari package datasets, dapat dibuat suatu model regresi dengan peubah respon Employed dan peubah penjelas semua kolom kecuali Employed. Hitunglah dugaan koefisien model regresi tersebut dengan

SWEEP Operator menggunakan R
Dekomposisi Cholesky menggunakan R
Menggunakan rumus Least-Squared

\[ b = (X'X)^{-1}X'y \]

Bandingkan waktu komputasi antara a,b dan c! Manakah yang lebih cepat?

Jawaban

Sebelum menjawab pertanyan-pertanyan nomor 3, kita import terlebih dahulu data nya

data("longley",package = "datasets")
dim(longley)

## [1] 16  7

longley

SWEEP Operator menggunakan R

Langkah Pertama: Buat matriks M

# membuat matrix X
X <- model.matrix(formula(Employed~.),data = longley)
dim(X)

## [1] 16  7

##      (Intercept) GNP.deflator     GNP Unemployed Armed.Forces Population Year
## 1947           1         83.0 234.289      235.6        159.0    107.608 1947
## 1948           1         88.5 259.426      232.5        145.6    108.632 1948
## 1949           1         88.2 258.054      368.2        161.6    109.773 1949
## 1950           1         89.5 284.599      335.1        165.0    110.929 1950
## 1951           1         96.2 328.975      209.9        309.9    112.075 1951
## 1952           1         98.1 346.999      193.2        359.4    113.270 1952
## 1953           1         99.0 365.385      187.0        354.7    115.094 1953
## 1954           1        100.0 363.112      357.8        335.0    116.219 1954
## 1955           1        101.2 397.469      290.4        304.8    117.388 1955
## 1956           1        104.6 419.180      282.2        285.7    118.734 1956
## 1957           1        108.4 442.769      293.6        279.8    120.445 1957
## 1958           1        110.8 444.546      468.1        263.7    121.950 1958
## 1959           1        112.6 482.704      381.3        255.2    123.366 1959
## 1960           1        114.2 502.601      393.1        251.4    125.368 1960
## 1961           1        115.7 518.173      480.6        257.2    127.852 1961
## 1962           1        116.9 554.894      400.7        282.7    130.081 1962
## attr(,"assign")
## [1] 0 1 2 3 4 5 6

# membuat vektor y
y <- as.matrix(longley$Employed)
dim(y)

## [1] 16  1

##         [,1]
##  [1,] 60.323
##  [2,] 61.122
##  [3,] 60.171
##  [4,] 61.187
##  [5,] 63.221
##  [6,] 63.639
##  [7,] 64.989
##  [8,] 63.761
##  [9,] 66.019
## [10,] 67.857
## [11,] 68.169
## [12,] 66.513
## [13,] 68.655
## [14,] 69.564
## [15,] 69.331
## [16,] 70.551

#membuat matriks M
M <- cbind(X,y)
dim(M)

## [1] 16  8

##      (Intercept) GNP.deflator     GNP Unemployed Armed.Forces Population Year
## 1947           1         83.0 234.289      235.6        159.0    107.608 1947
## 1948           1         88.5 259.426      232.5        145.6    108.632 1948
## 1949           1         88.2 258.054      368.2        161.6    109.773 1949
## 1950           1         89.5 284.599      335.1        165.0    110.929 1950
## 1951           1         96.2 328.975      209.9        309.9    112.075 1951
## 1952           1         98.1 346.999      193.2        359.4    113.270 1952
## 1953           1         99.0 365.385      187.0        354.7    115.094 1953
## 1954           1        100.0 363.112      357.8        335.0    116.219 1954
## 1955           1        101.2 397.469      290.4        304.8    117.388 1955
## 1956           1        104.6 419.180      282.2        285.7    118.734 1956
## 1957           1        108.4 442.769      293.6        279.8    120.445 1957
## 1958           1        110.8 444.546      468.1        263.7    121.950 1958
## 1959           1        112.6 482.704      381.3        255.2    123.366 1959
## 1960           1        114.2 502.601      393.1        251.4    125.368 1960
## 1961           1        115.7 518.173      480.6        257.2    127.852 1961
## 1962           1        116.9 554.894      400.7        282.7    130.081 1962
##            
## 1947 60.323
## 1948 61.122
## 1949 60.171
## 1950 61.187
## 1951 63.221
## 1952 63.639
## 1953 64.989
## 1954 63.761
## 1955 66.019
## 1956 67.857
## 1957 68.169
## 1958 66.513
## 1959 68.655
## 1960 69.564
## 1961 69.331
## 1962 70.551

Langkah Kedua: Hitung matriks A

A <- crossprod(M)
dim(A)

## [1] 8 8

##              (Intercept) GNP.deflator          GNP Unemployed Armed.Forces
## (Intercept)       16.000       1626.9     6203.175     5109.3       4170.7
## GNP.deflator    1626.900     167172.1   646700.650   528908.0     429317.4
## GNP             6203.175     646700.6  2553151.560  2065054.2    1663294.5
## Unemployed      5109.300     528908.0  2065054.181  1762542.7    1314528.0
## Armed.Forces    4170.700     429317.4  1663294.516  1314528.0    1159816.8
## Population      1878.784     192139.7   738680.235   606648.6     492386.4
## Year           31272.000    3180539.9 12131170.206  9990586.4    8153706.8
##                 1045.072     106816.2   410322.735   336197.8     274094.1
##               Population     Year            
## (Intercept)     1878.784    31272    1045.072
## GNP.deflator  192139.651  3180540  106816.177
## GNP           738680.235 12131170  410322.735
## Unemployed    606648.556  9990586  336197.802
## Armed.Forces  492386.424  8153707  274094.133
## Population    221340.143  3672577  123068.464
## Year         3672577.089 61121464 2042836.838
##               123068.464  2042837   68445.977

Langkah Kedua: Hitung koefisien Regresi dengan SWEEP operator

Karena banyaknya peubah penjelas ada 6 ditambah 1 intercept maka banyaknya sweep yang dilakukan untuk mendapatkan semua koefisien regresi adalah 7

res_sweep <- sweep.operator(A,k = 1:7)
dim(res_sweep)

## [1] 8 8

koefisien regresi bisa diambil dari kolom ke 8 baris 1 sampai 7. Hal ini sesuai dengan definisi matriks hasil akhir dari SWEEP Operator

\[ \text{SWEEP}(A,k) = \begin{pmatrix} (X'X)_{p\times p}^{-1} & |& b_{p\times 1} \\ -b_{1\times p}'& |& \text{SSE}_{1\times1} \end{pmatrix} \]

coef_sweep <- res_sweep[-8,8]
length(coef_sweep)

## [1] 7

coef_sweep

##   (Intercept)  GNP.deflator           GNP    Unemployed  Armed.Forces 
## -3.482259e+03  1.506187e-02 -3.581918e-02 -2.020230e-02 -1.033227e-02 
##    Population          Year 
## -5.110410e-02  1.829151e+00

Dekomposisi Cholesky menggunakan R

Langkah Pertama: Buat matriks X dan vektor y

# membuat matrix X
X <- model.matrix(formula(Employed~.),data = longley)
dim(X)

## [1] 16  7

# membuat vektor y
y <- as.matrix(longley$Employed)
dim(y)

## [1] 16  1

#membuat matriks M
M <- cbind(X,y)
dim(M)

## [1] 16  8

Langkah Kedua: Masukan X dan y tersebut ke fungsi lm_chol

coef_chol <- lm_chol(X,y)
length(coef_chol)

## [1] 7

coef_chol

## [1] -3.482259e+03  1.506187e-02 -3.581918e-02 -2.020230e-02 -1.033227e-02
## [6] -5.110411e-02  1.829151e+00

Menggunakan rumus Least-Squared

\[ b = (X'X)^{-1}X'y \]

Langkah Pertama : Buat matriks X dan vektor y

# membuat matrix X
X <- model.matrix(formula(Employed~.),data = longley)
dim(X)

## [1] 16  7

# membuat vektor y
y <- as.matrix(longley$Employed)
dim(y)

## [1] 16  1

#membuat matriks M
M <- cbind(X,y)
dim(M)

## [1] 16  8

Langkah Kedua : Hitung menggunakan rumus least squared

lm_lsq <- function(X,y){as.vector(solve(crossprod(X),crossprod(X,y)))}

coef_ls <- lm_lsq(X,y)
length(coef_ls)

## [1] 7

coef_ls

## [1] -3.482259e+03  1.506187e-02 -3.581918e-02 -2.020230e-02 -1.033227e-02
## [6] -5.110411e-02  1.829151e+00

Bandingkan waktu komputasi antara a sampai c! Manakah yang lebih cepat?

Kita akan membuat user-defined function untuk SWEEP operator agar perbandingan waktu komputasi lebih fair.

lm_sweep <- function(X,y){
  M <- cbind(X,y)
  A <- crossprod(M)
  num_col <- ncol(X)
  res_sweep <- sweep.operator(A,k = 1:num_col)
  return(res_sweep[-(num_col+1),(num_col+1)])
}

Untuk membandingkan kecepatan, kita akan menggunakan fungsi microbenchmark dari package microbenchmark

speed_test <- microbenchmark(
  lm_sweep(X,y),
  lm_chol(X,y),
  lm_lsq(X,y),
  times = 100 # 100 kali ulangan 
)

ggplot2::autoplot(speed_test)

## Coordinate system already present. Adding new coordinate system, which will
## replace the existing one.

Soal 4

Dengan menggunakan data US Crime, dapat dibuat suatu model regresi dengan peubah respon ViolentCrimesPerPop dan peubah penjelas semua kolom kecuali ViolentCrimesPerPop. Hitunglah dugaan koefisien model regresi tersebut dengan

SWEEP Operator menggunakan R
Dekomposisi Cholesky menggunakan R
QR-decomposition menggunakan fungsi qr.solve dalam R
Menggunakan rumus Least-Squared

\[ b = (X'X)^{-1}X'y \]

Bandingkan waktu komputasi antara a sampai d! Manakah yang lebih cepat?

Jawaban

Sebelum menjawab pertanyan-pertanyan nomor 4, kita import terlebih dahulu data nya

us_crime <- read.csv("data/us_crime_ed2.csv")
dim(us_crime)

## [1] 817 101

head(us_crime)

SWEEP Operator menggunakan R

Langkah Pertama : Buat matriks M

# membuat matrix X
X <- model.matrix(formula(ViolentCrimesPerPop~.),data = us_crime)
dim(X)

## [1] 817 101

# membuat vektor y
y <- as.matrix(us_crime$ViolentCrimesPerPop)
dim(y)

## [1] 817   1

#membuat matriks M
M <- cbind(X,y)
dim(M)

## [1] 817 102

Langkah Kedua : Hitung matriks A

A <- crossprod(M)
dim(A)

## [1] 102 102

Langkah Kedua : Hitung koefisien Regresi dengan SWEEP operator

Karena banyaknya peubah penjelas ada 100 ditambah 1 intercept maka banyaknya sweep yang dilakukan untuk mendapatkan semua koefisien regresi adalah 101

res_sweep <- sweep.operator(A,k = 1:101)
dim(res_sweep)

## [1] 102 102

koefisien regresi bisa diambil dari kolom ke 102 baris 1 sampai 101. Hal ini sesuai dengan definisi matriks hasil akhir dari SWEEP Operator

\[ \text{SWEEP}(A,k) = \begin{pmatrix} (X'X)_{p\times p}^{-1} & |& b_{p\times 1} \\ -b_{1\times p}'& |& \text{SSE}_{1\times1} \end{pmatrix} \]

coef_sweep <- res_sweep[-102,102]
length(coef_sweep)

## [1] 101

head(coef_sweep,10)

##   (Intercept)    population householdsize  racepctblack  racePctWhite 
##   -0.05748964    1.15370365   -0.03192869   -0.15581456   -0.41394505 
##  racePctAsian   racePctHisp   agePct12t21   agePct12t29   agePct16t24 
##   -0.05582214   -0.20301467    0.31691329    0.03486022   -0.42825586

Dekomposisi Cholesky menggunakan R

Langkah Pertama : Buat matriks X dan vektor y

# membuat matrix X
X <- model.matrix(formula(ViolentCrimesPerPop~.),data = us_crime)
dim(X)

## [1] 817 101

# membuat vektor y
y <- as.matrix(us_crime$ViolentCrimesPerPop)
dim(y)

## [1] 817   1

Langkah Kedua : Masukan X dan y tersebut ke fungsi lm_chol

coef_chol <- lm_chol(X,y)
length(coef_chol)

## [1] 101

head(coef_chol,10)

##  [1] -0.05748964  1.15370365 -0.03192869 -0.15581456 -0.41394505 -0.05582214
##  [7] -0.20301467  0.31691329  0.03486022 -0.42825586

QR-decomposition menggunakan fungsi qr.solve dalam R

Langkah Pertama : Buat matriks X dan vektor y

# membuat matrix X
X <- model.matrix(formula(ViolentCrimesPerPop~.),data = us_crime)
dim(X)

## [1] 817 101

# membuat vektor y
y <- as.matrix(us_crime$ViolentCrimesPerPop)
dim(y)

## [1] 817   1

Langkah Kedua : Masukan X dan y tersebut ke fungsi qr.solve

coef_qr <- qr.solve(X,y)
length(coef_qr)

## [1] 101

head(as.vector(coef_qr),10)

##  [1] -0.05748964  1.15370365 -0.03192869 -0.15581456 -0.41394505 -0.05582214
##  [7] -0.20301467  0.31691329  0.03486022 -0.42825586

Menggunakan rumus Least-Squared

\[ b = (X'X)^{-1}X'y \]

Langkah Pertama : Buat matriks X dan vektor y

# membuat matrix X
X <- model.matrix(formula(ViolentCrimesPerPop~.),data = us_crime)
dim(X)

## [1] 817 101

# membuat vektor y
y <- as.matrix(us_crime$ViolentCrimesPerPop)
dim(y)

## [1] 817   1

Langkah Kedua : Hitung menggunakan rumus least squared

lm_lsq <- function(X,y){as.vector(solve(crossprod(X),crossprod(X,y)))}

coef_ls <- lm_lsq(X,y)
length(coef_ls)

## [1] 101

head(coef_ls,10)

##  [1] -0.05748964  1.15370365 -0.03192869 -0.15581456 -0.41394505 -0.05582214
##  [7] -0.20301467  0.31691329  0.03486022 -0.42825586

Bandingkan waktu komputasi antara a sampai d! Manakah yang lebih cepat?

Membuat user-defined function untuk SWEEP operator agar perbandingan waktu komputasi lebih fair.

lm_sweep <- function(X,y){
  M <- cbind(X,y)
  A <- crossprod(M)
  num_col <- ncol(X)
  res_sweep <- sweep.operator(A,k = 1:num_col)
  return(res_sweep[-(num_col+1),(num_col+1)])
}

Untuk membandingkan kecepatan, kita akan menggunakan fungsi microbenchmark dari package microbenchmark

speed_test <- microbenchmark(
  lm_sweep(X,y),
  lm_chol(X,y),
  qr.solve(X,y),
  lm_lsq(X,y),
  times = 100 # 100 kali ulangan 
)

summary(speed_test)

ggplot2::autoplot(speed_test)

## Coordinate system already present. Adding new coordinate system, which will
## replace the existing one.

Referensi

https://www.theissaclee.com/post/linearqrandchol/
https://blogs.sas.com/content/iml/2018/04/18/sweep-operator-sas.html
Goodnight, J. (1979). “A Tutorial on the SWEEP Operator.” The American Statistician.

Statistika dan Sains Data, IPB University, alfanugraha@apps.ipb.ac.id ↩︎