Analisis regresi merupakan salah satu bidang ilmu statistika yang umumnya digunakan untuk menyelidiki pola hubungan fungsional antara variabel respon dan variabel prediktor. Terdapat tiga model pendekatan regresi, yaitu regresi parametrik, regresi nonparametrik, dan regresi semiparametrik. Salah satu bentuk model regresi nonparametrik adalah model aditif jika semua fungsi pada prediktor aditif merupakan fungsi mulus. Selain itu, model aditif juga dapat dibentuk menjadi model regresi semiparametrik. Penggunaan model aditif hanya terbatas pada data yang berdistribusi normal, sehingga dikembangkan model aditif tergeneralisir (Generalized Additive Models atau GAM) untuk menangani kondisi tersebut, karenanya model aditif tergeneralisir lebih fleksibel dibanding model aditif. Variabel respon pada model aditif tergeneralisir diasumsikan berdistribusi keluarga eksponensial dan hubungannya dengan variabel prediktor tidak harus linier. Pencocokan model aditif tergeneralisir menggunakan metode penghalusan berupa penghalus diagram pencar untuk mengestimasi fungsi mulus, dan algoritma P-IRLS (Penalized Iteratively Re-weighted Least Square) untuk mengestimasi model. Adapun metode penghalusan yang dapat memberikan suatu hasil analisis numerik yang lebih baik menurut Hastie (1990) adalah penghalus spline. Penghalus spline mempunyai banyak kelebihan dibandingkan penghalus lainnya, salah satunya yaitu dapat mengatasi pola data yang menunjukkan naik/turun yang tajam serta kurva yang dihasilkan relatif mulus. Oleh karenanya, banyak peneliti yang menggunakan model analisis tergeneralisir dengan penghalus spline untuk menganalisis data. Dalam menganalisis model aditif tergeneralisir dengan penghalus spline, Wood (2006) telah memperkenalkan paket mgcv pada program R yang sebagian besar menggunakan pendekatan skrip. Hal ini tidak mudah bagi peneliti yang bukan berlatar belakang statistika maupun program statistika. Sehingga, dalam program ini paket mgcv pada program R untuk menganalisis model aditif tergeneralisir dengan penghalus spline dibuat dengan sistem online. Dengan adanya program ini dapat memudahkan peneliti yang kurang mengerti program statistika dalam melakukan analisis data, yaitu hanya dengan memasukkan data dan memilih menu secara online tanpa harus menginstal program R telebih dahulu.
Pada akhir kegiatan, pengguna diharapkan:
Analisis regresi merupakan salah satu bidang statistika yang memainkan peran sangat penting, umumnya digunakan untuk menyelidiki pola hubungan fungsional anatara variabel respon (variabel dependen) dan variabel prediktor (variabel independen). Dewasa ini, terdapat tiga model pendekatan regresi yang banyak dikembangkan oleh para peneliti, yaitu regresi parametrik, regresi nonparametrik, dan regresi semiparametrik.
dengan $x_i$ adalah variabel prediktor; $y_i$ adalah variabel respon; $f$ adalah kurva regresi; dan $\varepsilon_i$ adalah galat. Pada regresi parametrik terdapat asumsi yang sangat kaku dan kuat yaitu bentuk kurva regresinya diketahui (linier, kuadratik, kubik, polinomial derajat-p, eksponensial, dan lain-lain).
Berbeda dengan regresi parametrik yang memiliki asumsi yang sangat kaku dan kuat, regresi nonparametrik sangatlah fleksibel dan sangat obyektif. Pada regresi nonparametrik, biarkan data sendiri yang akan mencari bentuk estimasi dari kurva regresinya, tanpa dipengaruhi oleh faktor subyektifitas si peneliti seperti pada regresi parametrik yang cenderung ada unsur pemaksaan dari peneliti. Dalam regresi nonparametrik, bentuk kurva regresinya tidak diketahui dan hanya diasumsikan mulus (smooth). Model regresi nonparametrik dapat dituliskan seperti berikut
$$\begin{align*}y_i=s(t_i)+\varepsilon_i, i=1,2,...,n\end{align*}$$dengan $t_i$ adalah variabel prediktor; $y_i$ adalah variabel respon; $s$ adalah kurva regresi dan $\varepsilon_i$ adalah galat.
Selain model pendekatan regresi parametrik dan regresi nonparametrik, terdapat pula model pendekatan regresi semiparametrik yang merupakan gabungan keduanya. Para statistikawan memandang kurva regresi dapat diklasifikasikan ke dalam dua komponen, yaitu komponen parametrik (bentuk kurva regresinya diketahui) dan komponen nonparametrik (bentuk kurva regresinya tidak diketahui). Regresi semiparametrik ini digunakan untuk menangani data yang sebagian datanya memiliki bentuk kurva regresi yang diketahui dan sebagian lainnya tidak diketahui bentuk kurvanya. Model regresi semiparametrik secara umum dapat dituliskan seperti berikut
$$\begin{align*}y_i=f(x_i)+s(t_i)+\varepsilon_i,i=1,2,...,n\end{align*}$$dengan $f$ adalah kurva regresi parametrik dan $s$ adalah kurva regresi nonparametrik.
Model aditif tergeneralisir atau biasa disebut dengan GAM (Generalized Additive Models) ini mengembangkan model aditif $Y=\alpha+\sum_{j=1}^{p}f_{j}(X_{j})+\varepsilon$, ke dalam bentuk distribusi keluarga eksponensial. Model aditif tergeneralisir juga dapat disebut sebagai model pengembangan dari model linier tergeneralisir, karena model ini mengembangkan model linier tergeneralisir dengan cara yang sama dengan model aditif mengembangkan model linier, yaitu mengganti prediktor linier $\sum_{j=1}^{p}\beta_{j}X_{j}$ dengan prediktor aditif $\sum_{j=1}^{p}f_{j}(X_{j})$. Hal tersebut membuat model aditif tergeneralisir menjadi lebih fleksibel daripada model linier tergeneralisir maupun model aditif.
Seperti pada model linier tergeneralisir, model aditif tergeneralisir juga memuat komponen acak, komponen tetap yaitu prediktor aditif, dan fungsi link yang menghubungkan dua komponen tersebut. Variabel respon $Y$ sebagai komponen acak diasumsikan memiliki fungsi kepadatan keluarga eksponensial sebagai berikut:
$$\begin{align*}f_Y(y;\theta;\phi)=\exp\left[\frac{y\theta-b(\theta)}{a(\phi)}+c(y,\phi)\right]\end{align*}$$dengan $\theta$ adalah parameter natural dan $\phi$ adalah parameter skala. Hubungan antara mean $\mu$ dari variabel respon dan komponen aditif $\eta(X)=s_0+\sum_{j=1}^{p}s_j(X_j)$ didefinisikan oleh fungsi link $g(\mu)=\eta$ dan secara umum fungsi link yang paling umum digunakan adalah link kanonik dengan $\eta=\theta$ (Hastie dan Tibshirani, 1990). Model aditif tergeneralisir dapat diformulasikan sebagai berikut:
$$\begin{align*}g(\mu)=s_0+\sum_{j=1}^{p}s_j(X_j)+\varepsilon\end{align*}$$dengan variabel responnya diasumsikan berdistribusi keluarga eksponensial dan $s(.)$ merupakan fungsi mulus.
Model aditif tergeneralisir dan model linier tergeneralisir dapat diterapkan dalam situasi yang mirip, tetapi kedua model tersebut memberikan tujuan analitik yang berbeda. Model linier tergeneralisir menekankan estimasi dan inferensi untuk parameter model, sementara model aditif tergeneralisir fokus pada eksplorasi data nonparametrik (Statsoft, 2016).
Pada model aditif tergeneralisir, fungsi atau komponen nonparametrik diestimasi oleh fungsi penghalus yakni penghalus diagram pencar (scatterplot smoother). Hastie dan Tibshirani (1990), menjelaskan berbagai macam penghalus diagram pencar, salah satunya adalah penghalus spline. Misalkan diberikan data ${(x_i,y_i ),i=1,\cdots,n}$, dengan $n$ adalah banyaknya pengamatan dan hubungan antara $x_i$ dan $y_i$ diasumsikan mengikuti model regresi nonparametrik berikut
$$y_i=f(x_i)+\varepsilon_i,x_i\in[a,b]$$dimana dalam hal ini $f(x_i)$ adalah fungsi yang tidak diketahui bentuknya dan hanya diasumsikan mulus. Menurut Herawati (2011), penghalus spline akan mengestimasi fungsi $f$ sebagai solusi dari masalah optimasi yaitu dengan mencari $\hat{f}\in L_2[a,b]$ yang meminimumkan persamaan penalized least squares berikut
$$S(f)=\sum_{i=1}^{n}(y_i-f(x_i))^2+\lambda\int_a^b(f^{''}(x))^2dx$$dimana $L_2[a,b]$ menyatakan himpunan fungsi kuadrat terintegral pada interval $[a,b]$; $\lambda$ adalah nilai konstan (fixed constant); dan $a\le x_i\le\cdots\le x_n\le b$. Suku pertama pada persamaan $(2)$ digunakan untuk mengukur kerapatan data, yang tidak lain merupakan jumlah kuadrat sisaan, sedangkan suku kedua merupakan penalti kekasaran (roughness penalty) yang memberikan ukuran kemulusan atau kekasaran kurva dalam memetakan data.Parameter $\lambda$ merupakan parameter penghalus (smoothing parameter) yang bernilai positif. Nilai $\lambda$ yang besar akan menghasilkan kurva yang mulus, sedangkan nilai $\lambda$ yang kecil akan menghasilkan kurva yang kasar. Parameter penghalus ini memainkan peran sentral dalam mengendalikan keseimbangan (trade off) antara ketepatan model (goodness of fit) dan kemulusan kurva (Herawati, 2011).
Salah satu penghalus spline adalah spline kubik dengan titik perubahan yang terjadi di dalam suatu kurva yang disebut titik knot. Fungsi spline berorde $m$ dengan titik-titik knot $K_1,K_2,\cdots,K_N$ secara umum dapat disajikan dalam bentuk (Herawati, 2011):
Jika pada persamaan $(3)$ diambil nilai $m=4$, maka diperoleh spline kubik. Menurut Green dan Silverman dalam Herawati (2011), suatu fungsi $f(x)$ disebut penghalus spline kubik jika memenuhi dua kondisi, yaitu kondisi pertama, fungsi $f(x)$ pada setiap interval $(a,x_1),(x_1,x_2),(x_2,x_3),\cdots,(x_n,b)$ adalah polinomial kubik, dan kondisi kedua yaitu fungsi $f(x)$ kontinu pada turunan pertama dan kedua pada setiap $x_i\in[a,b]$ dimana $x_i$ merupakan titik knot.
Misalkan terdapat sebuah model yang mengandung satu fungsi mulus seperti pada model $(1)$ di atas. Maka estimasi fungsi mulus $f$ dengan regresi spline kubik dapat dinyatakan ke dalam bentuk
$$f(x)=\sum_{j=1}^k h_j(x)\beta_j$$dimana $k$ adalah banyaknya knot; $h_j(x)$ merupakan fungsi basis ke-$j$ untuk regresi spline kubik; dan $\beta_j$ adalah parameter yang tidak diketahui (Wood, 2006). Menurut Wood (2006), untuk sebarang fungsi $f$ yang mempunyai sebuah ekspansi basis seperti pada $(4)$, memungkinkan untuk menulis penalti kekasaran seperti berikut
$$\int_{x_1}^{x_k}(f^{''}(x))^2 dx=\boldsymbol{\beta}^T \boldsymbol{S\beta}$$dengan $x_1,\cdots,x_k$ adalah titik knot; $\boldsymbol{\beta}^T=(\beta_1,\cdots,\beta_k)$; dan $S$ adalah matriks koefisien yang diketahui. Matriks penalti untuk basis regresi spline kubik ini adalah
$$\begin{align*}\boldsymbol{S}=\boldsymbol{D}^T \boldsymbol{B}^{-1} \boldsymbol{D}\end{align*}$$dengan $\boldsymbol D$ adalah matriks upper triagonal dan $\boldsymbol B$ adalah matriks symmetric tridiagonal yang didefinisikan pada tabel berikut.
| $D_{j,j}=1/h_j$ | $D_{j,j+1}=-1/h_j-1/h_{j+1}$ | $D_{j,j+2}=1/h_{j+1}$ |
| $B_{j,j}=(h_j+h_{j+1})/3$ | $ $ | $j=1\cdots k-2$ |
| $ $ | $ $ | $ $ |
| $B_{j,j+1}=h_{j+1}/6$ | $B_{j+1,j}=h_{j+1}/6$ | $j=1\cdots k-3$ |
Fungsi regresi spline kubik pada $(4)$ dapat dituliskan ke dalam bentuk matriks yaitu
$$\boldsymbol{f}=\boldsymbol{H\beta}$$dimana $\boldsymbol{f}=(f(x_1),\cdots,f(x_n))^T$; $\boldsymbol{\beta}=(\beta_1,\cdots,\beta_n)^T$; dan $H_{i,j}=h_j(x_i)$. Nilai estimasi parameter $\beta$ pada fungsi mulus $f$ diperoleh dengan meminimumkan persamaan $(2)$. Persamaan $(2)$ di atas dapat dituliskan dalam bentuk matriks seperti berikut
$$S(f)=(\boldsymbol{y-H\beta})^T (\boldsymbol{y-H\beta})+ \lambda\boldsymbol{\beta}^T \boldsymbol{S\beta}$$Syarat perlu agar persamaan di atas minimum adalah $\frac{\partial {S(f)}}{\partial \beta}=0$, sehingga diperoleh nilai estimasi untuk parameter $\beta$ yaitu
$$\boldsymbol{\beta}=(\boldsymbol{H}^T \boldsymbol{H}+\lambda\boldsymbol{S})^{-1} \boldsymbol{H}^T \boldsymbol{y}$$dengan mensubtitusikan persamaan $(8)$ ke $(6)$, maka didapatkan estimasi fungsi mulus $f$ dengan penghalus regresi spline kubik yaitu
$$\boldsymbol{f}=\boldsymbol{H}(\boldsymbol{H}^T\boldsymbol{H}+\lambda \boldsymbol{S})^{-1} \boldsymbol{H}^T \boldsymbol{y}$$dimana $\boldsymbol{S}=\boldsymbol{D}^T \boldsymbol{B}^{-1} \boldsymbol{D}$.
Regresi spline kubik siklik adalah spesifikasi dari regresi spline kubik yang nilai fungsi dan turunannya hingga order kedua di dua titik ujungnya adalah sama. Seperti regresi spline kubik yang telah dijelaskan sebelumnya, estimasi fungsi mulus $f$ pada model $(1)$ dengan regresi spline kubik siklik juga dapat dinyatakan ke dalam bentuk berikut
$$f(x)=\sum_{j=1}^{k-1}g_j(x)\beta_j$$dimana $k$ adalah banyaknya knot; $g_j(x)$ merupakan fungsi basis ke-$j$ untuk regresi spline kubik siklik; dan $\beta_j$ adalah parameter yang tidak diketahui (Wood, 2006). Berbeda dengan regresi spline kubik sebelumnya, pada regresi spline kubik siklik $\beta_1=\beta_k$ dan $g_1(x)=g_k(x)$, karenanya $j$ pada persamaan $(10)$ bergerak dari $1$ hingga $k-1$. Fungsi $f$ pada $(1)$ mempunyai sebuah ekspansi basis, sehingga memungkinkan untuk menulis penalti kekasaran (roughness penalty) seperti pada persamaan $(5)$ berikut
$$\begin{align*}\int_{x_1}^{x_k}(f^{''}(x))^2 dx=\boldsymbol{\beta}^T \boldsymbol{S\beta}\end{align*}$$dalam hal ini $\boldsymbol{\beta}^T=(\beta_1,\cdots,\beta_{k-1})$; dan $S$ adalah matriks koefisien yang diketahui. Matriks penalti untuk basis regresi spline kubik siklik ini adalah
$$\begin{align*}\boldsymbol{S}=\boldsymbol{M}^T \boldsymbol{N}^{-1} \boldsymbol{M}\end{align*}$$dengan ${\boldsymbol{M}}$ dan ${\boldsymbol{N}}$ adalah matriks yang didefinisikan pada tabel berikut.
| $N_{j-1,j}=N_{j,j-1}=h_{j-1}/6$ | $N_{j,j}=(h_{j-1}+h_1)/3$ | $ $ |
| $M_{j-1,j}=M_{j,j-1}=1/h_{j+1}$ | $M_{j,j}=-1/h_{j-1}-1/h_j$ | $j=2\cdots k-1$ |
| $ $ | $ $ | $ $ |
| $N_{1,1}=(h_{k-1}+h_1)/3$ | $N_{1,k-1}=h_{k-1}/6$ | $N_{k-1,1}=h_{k_1}/6$ |
| $M_{1,1}=-1/h_1-1/h_{k-1}$ | $M_{1,k-1}=1/h_{k-1}$ | $M_{k-1,1}=1/h_{k_1}$ |
Seperti pada regresi spline kubik sebelumnya, fungsi regresi spline kubik siklik pada $(10)$ juga dapat dituliskan ke dalam bentuk matriks, yaitu
$$\boldsymbol{f}=\boldsymbol{G\beta}$$dimana $\boldsymbol{f}=(f(x_1),\cdots,f(x_n))^T$; $\boldsymbol{\beta}=\beta(\beta_1,\cdots,\beta_{k-1})^T$; dan $G_{i,j}=g_j (x_i)$. Dengan cara yang sama dapat diperoleh nilai estimasi untuk parameter $\beta$ yaitu
$$\boldsymbol{\beta}=(\boldsymbol{G}^T \boldsymbol{G}+\lambda \boldsymbol{S})^{-1} \boldsymbol{G}^T \boldsymbol{y}$$dengan mensubtitusikan persamaan $(12)$ ke $(11)$, maka didapatkan estimasi fungsi mulus $f$ dengan penghalus regresi spline kubik siklik yaitu
$$\boldsymbol{f}=\boldsymbol{G}(\boldsymbol{G}^T \boldsymbol{G}+\lambda \boldsymbol{S})^{-1} \boldsymbol{G}^T \boldsymbol{y}$$dimana $\boldsymbol{S}=\boldsymbol{M}^T \boldsymbol{N}^{-1} \boldsymbol{M}$.
Berdasarkan sifat yang dimiliki di atas, maka regresi spline kubik siklik dapat digunakan untuk menganalisis data time series yang berpola periodik atau siklik.
Menurut Budiantara (dalam Tripena, 2011), bentuk estimator spline sangat dipengaruhi oleh nilai parameter penghalus $\lambda$. Pemilihan $\lambda$ optimal dalam regresi spline pada hakikatnya merupakan pemilihan lokasi titik-titik knot. Mengingat tujuan dari pendekatan regresi nonparametrik, yaitu ingin didapatkan kurva mulus yang mempunyai $\lambda$ optimal menggunakan data amatan sebanyak $n$, maka diperlukan ukuran kinerja atas estimator yang dapat diterima secara universal. Eubank (1999) menyebutkan bahwa, ukuran kinerja atas estimator tersebut adalah Generalized Cross-Validation (GCV) yang diformulasikan sebagai berikut
$$GCV(\lambda)=\frac{1}{n} \sum_{i=1}^n\left(\frac{y_i-\hat{f}_{\lambda}(x_i)} {1-tr({S}_{ii}(\lambda))}\right)^2$$dengan $\boldsymbol{\hat{f}_{\lambda}}=\boldsymbol{S}_{\lambda}+ \boldsymbol{y}$ dan $S_{ii}(\lambda)$ adalah elemen diagonal dari $\boldsymbol{S}_{\lambda}$, dimana $\boldsymbol{S}_{\lambda}$ adalah matriks penghalus yang berukuran $n \times n$. Kriteria GCV diharapkan memiliki nilai yang minimum, sehingga model regresi spline dapat dikatakan memiliki nilai $\lambda$ yang optimal.
Selain nilai GCV, nilai statistik AIC (Akaike's Information Criterion) juga diperlukan untuk memilih ketepatan atau kecocokan model secara keseluruhan. Besarnya AIC dapat dilihat pada persamaan berikut
$$AIC=-2l(\boldsymbol{\hat{\theta}})+2q$$dengan $l(\boldsymbol{\hat{\theta}})$ adalah nilai likelihood dari model dan $q$ adalah banyaknya parameter dalam model. Nilai AIC yang lebih kecil menunjukkan model yang lebih tepat atau cocok (Tirta, 2009).
Adapun langkah-langkah kerja analisis data menggunakan GAM dengan penghalus spline pada kegiatan ini adalah sebagai berikut:
|
Khusus untuk Impor Data, Pilih File: |
Terdapat dua cara yang dapat digunakan untuk memilih komponen parametrik dan nonparametrik, yaitu dengan melihat pola hubungan antar variabel pada scatterplot dan melihat nilai koefisien determinasi ($R^2$).
1. ScatterplotScatterplot matrix dapat memberikan deskripsi mengenai hubungan kelinieran antar variabel. Pada scatterplot matrix, variabel respon dibaca secara vertikal dan variabel prediktor dibaca secara horisontal. Terdapat dua garis dengan warna yang berbeda pada scatterplot, yaitu hijau dan merah. Garis hijau merupakan garis linier, sedangkan garis merah merupakan garis regresi yang dihasilkan dari plot data yang ada. Ketika variabel respon dan variabel prediktor diinteraksikan dan menghasilkan garis regresi yang mendekati garis linier yang terbentuk, maka variabel prediktor tersebut termasuk dalam komponen paramterik. Sedangkan apabila garis regresi yang terbentuk menjauhi garis linier, maka variabel prediktor tersebut termasuk dalam komponen nonparametrik.
Variabel:Catatan: Pilih minimal 2 (dua) variabel numerik dengan cara menekan tombol ctrl dan klik variabel-variabel yang dipilih.
|
| |
Koefisien determinasi dalam regresi linier sering diartikan sebagai seberapa besar kemampuan variabel prediktor dalam menjelaskan varians dari variabel responnya. Koefisien determinasi juga dapat menjadi indikator kebaikan suatu model. Dalam hal ini, nilai koefisien determinasi didapat dengan memodelkan variabel respon ($Y$) dengan variabel prediktor ($X$) menggunakan model linier. Sehingga, nilai $R^2$ yang lebih besar menunjukkan variabel prediktor yang lebih cocok jika digunakan pendekatan linier yang kemudian dimodelkan secara parametrik, dan begitu pula sebaliknya.
Pilih Variabel respon ($Y$) dan $X$ untuk praanalisis (menentukan kandidat nonparametrik)| Y: | X: |
Distribusi yang digunakan dalam model aditif tergeneralisir adalah distribusi keluarga eksponensial. Terdapat 4 (empat) distribusi yang disediakan dalam program ini, yaitu distribusi normal, gamma, poisson, dan binomial. Untuk data kontinu, dapat menggunakan distribusi normal atau gamma. Sedangkan untuk data diskrit (cacahan), dapat menggunakan distribusi poisson atau binomial. Pemilihan distribusi yang cocok, dapat dilihat dari nilai AIC yang minimum dan histogram distribusinya. Pada histogram distribusi, distribusi dikatakan cocok jika garis merah yang menggambarkan probability distribution function (pdf) mengikuti histogram frekuensi dari data variabel respon ($Y$).
Pilih Distribusi (tanpa link)
|
|
Selain output di atas, terdapat beberapa output lain yang dapat dipilih mengenai hasil pencocokan model menggunakan GAM seperti berikut.
Pilih Output:Uji signifikansi dilakukan untuk mengetahui variabel-variabel prediktor yang signifikan atau memberikan pengaruh terhadap variabel respon. Uji signifikansi ini dapat dilakukan dengan menggunakan nilai p-value. Variabel prediktor dikatakan signifikan jika nilai p-value$\leq\alpha$, dan begitupun sebaliknya. Dimana $\alpha$ adalah taraf signifikansi yang bernilai 0,01, 0,05, atau 0,1. Nilai p-value ini diperoleh dan dapat dilihat pada bagian hasil pencocokan model di atas.
Model aditif tergeneralisir (GAM) memiliki asumsi yang lebih fleksibel jika dibandingkan dengan model linier maupun model linier tergeneralisir (GLM), yaitu distribusi yang digunakan merupakan distribusi keluarga eksponensial dan hubungan antara variabel respon dengan variabel prediktor tidak harus linier. Fungsi mulus pada GAM diestimasi menggunakan penghalus diagram pencar, salah satu penghalus yang sering digunakan para peneliti karena kelebihan yang dimilikinya adalah penghalus spline. Dalam GAM ini, model terbaik dapat diketahui melalui nilai AIC yang minimum. Selain itu, model dengan parameter penghalus optimal dapat dilihat melalui nilai GCV yang minimum.
Eubank, W. 1999. Nonparametric Regression and Spline Smoothing. New York: Marcel Dekker.
Hastie, T. J., dan Tibshirani, R. J. 1990. Generalized Additive Models. London: Chapman & Hall.
Statsoft. 2016. Generalized Additive Models. http://www.statsoft.com/textbook/generalized-additive-models [23 Maret 2016].
Tirta, I. M. 2009. Analisis Regresi dengan R. Jember: Universitas Jember.
Wood, S. N. 2006. Generalized Additive Models: An Introduction with R. London: Chapman & Hall/CRC Press.