Machine Learning: Multiple Linear Regression

Catatan penting : Jika Anda benar-benar awam tentang apa itu Python, silakan klik artikel saya ini. Jika Anda awam tentang R, silakan klik artikel ini.

Sebelumnya kita sudah bersama-sama belajar tentang simple linear regression (SLR), kali ini kita belajar yang sedikit lebih advanced yaitu multiple linear regression (MLR). Apa perbedaannya? Pada pembahasan SLR, kita memiliki satu variabel dependen (gaji yang diterima per tahun), dan satu variabel independen (lamanya bekerja). Jika kita ingin mengetahui pengaruh variabel dependen (output) terhadap lebih dari satu variabel independen (data input, atau disebut juga dengan prediktor), maka kita harus menggunakan MLR.

Persamaan SLRy = a_{{0}} + a_{{1}}*x_{{1}}

di mana y adalah variabel dependen, a_{{0}} adalah konstanta, a_{{1}} adalah koefisien untuk x_{{1}}, dan x_{{1}} adalah variabel independen

Persamaan MLRy = a_{{0}} + a_{{1}}*x_{{1}} + a_{{2}}*x_{{2}} + ... + a_{{n}}*x_{{n}}

di mana MLR memiliki fungsi yang mirip dengan SLR namun memiliki variabel independen sebanyak n

Misal, dalam pembelajaran kali ini ada dataset keuntungan yang dimiliki oleh 50 perusahaan baru yang berada di US. Kita ingin melihat apakah ada hubungan keuntungan yang didapat (profit) terhadap pengeluaran-pengeluaran dan wilayah domisili perusahaan-perusahaan ini. Ada 4 variabel independen yang kita miliki yaitu data pengeluaran R & D (riset dan development), data pengeluaran marketing, data pengeluaran administrasi (seperti gaji, tagihan, dan lain-lain), dan wilayah di mana 50 perusahaan ini berada. Karena ada lebih dari 1 variabel independen, maka kita harus menggunakan multiple linear regression. Kita juga penasaran, dari 4 variabel independen ini mana yang memiliki kontribusi terbesar terhadap profit. Jika kita sudah mengetahui variabel independen mana yang memiliki kontribusi terbesar terhadap profit, maka kita juga bisa mengambil keputusan untuk memprioritaskan investasi dan fokus kita di variabel-variabel independen ini. Tentunya masih banyak lagi keputusan-keputusan yang bisa diambil dari model regresi kita kali ini.

Perlu diperhatikan, bahwa saat kita menggunakan MLR, kita dari awal sudah mengasumsikan bahwa hubungan keempat variabel independen terhadap variabel dependen adalah linear. Tentu saja, kenyataannya hubungan yang kita temukan setelah mengevaluasi model kita nanti bisa saja tidak linear. Dan jika hal ini terjadi, maka model MLR tidak cocok untuk mencari hubungan ketiganya. Tenang saja, untuk contoh kali ini, datanya sudah disesuaikan, sehingga ada hubungan linear antara variabel dependen dan independen.

Hal lain yang perlu diperhatikan adalah dalam membuat model linear regresi (simple atau multi), maka variabel dependen harus memenuhi beberapa asumsi lainnya. Di mana jika asumsi ini tidak terpenuhi maka model kita dinyatakan tidak cukup baik. Asumsi-asumsi tersebut antara lain:

  • Linearity
  • Homoscedasticity
  • Multivariate normality
  • Independence of errors
  • Lack of multicollinearity

Walau demikian, kita tidak akan membahas asumsi-asumsi di atas (akan di bahas di topik lainnya). Untuk mempermudah, kita asumsikan bahwa semua asumsi tersebut terpenuhi.

Baik, mari kita mulai membangun model machine learning MLR kita. Langkah pertama, Anda harus mendownload datasetnya di link ini.

Mari kita mulai dengan mendata variabel apa saja yang kita miliki.

  • Variabel dependen = Profit (data numerik)
  • Variabel independen 1 = Biaya R&D (data numerik)
  • Variabel independen 2 = Biaya administrasi (data numerik)
  • Variabel independen 3 = Biaya marketing (data numerik)
  • Variabel independen 4 = Wilayah perusahaan (data kategori)

Semua satuan untuk variabel independen 1 sampai 4 dan variabel dependen adalah dalam dollar ($) per tahun.

 

 

Pada jenis data variabel yang dimiliki, kita melihat bahwa variabel independen 4 tidak memiliki jenis data numerik. Padahal model regresi (simple atau multi) adalah model di mana semua datanya haruslah numerik. Lalu apa yang harus dilakukan?

Untuk mensiasati ini, maka kita perlu membuat dummy variable. Dummy variable adalah konversi data kategori ke bentuk numerik yang terdiri dari dua nilai yaitu 0 dan 1. Ilustrasi dari dummy variabel ini bisa dilihat pada tabel di bawah ini:

 

 

Bisa Anda lihat, bahwa kita merubah data kategori menjadi boolean (True or False). Nilainya bernilai True (1) jika datanya sama. Misal, baris pertama New York bernilai satu untuk kolom New York, namun FALSE (0) untuk California dan Florida. Begitu seterusnya sampai baris-baris lainnya. Namun perlu diingat, jika kita memiliki 3 kategori seperti contoh di atas, maka total dummy variabel yang dibutuhkan adalah n-1, dalam hal ini 3-1 = 2 dummy variabel. Hal ini dilakukan untuk menghindari jebakan dummy variabel (dummy variable trap).

Kira-kira funsgi multiple regresinya akan tampak sebagai berikut:

 

Profit = a_{{0}} + a_{{1}}*R.D + a_{{2}}*Marketing + a_{{3}}*Administrasi + a_{{4}}*Wilayah1 + a_{{5}}*Wilayah2

 

Tips: Saat membuang salah satu dummy variabel, tidak ada aturan apakah harus dummy variabel pertama, kedua atau ketiga, sehingga pilihan bebas terserah Anda. Walau demikian, pada contoh formula di atas, kita membuang dummy variabel ketiga (Wilayah3).

Mari sekarang kita buat modelnya step-by-step dilengkapi penjelasannya.

 

 

Bahasa Python

# Mengimpor library yang diperlukan
import numpy as np
import matplotlib.pyplot as plt
import pandas as pd

# Mengimpor dataset
dataset = pd.read_csv('50_Startups.csv')
X = dataset.iloc[:, :-1].values
Tampilkan_X = pd.DataFrame(X)
y = dataset.iloc[:, 4].values

# Encode categorical data
from sklearn.preprocessing import LabelEncoder, OneHotEncoder
labelencoder_X = LabelEncoder()
X[:, 3] = labelencoder.fit_transform(X[:, 3])
onehotencoder = OneHotEncoder(categorical_features = [3])
X = onehotencoder.fit_transform(X).toarray()

# Menghindari jebakan dummy variabel
X = X[:, 1:]

# Membagi data menjadi the Training set and Test set
from sklearn.cross_validation import train_test_split
X_train, X_test, y_train, y_test = train_test_split(X, y, test_size = 0.2, random_state = 0)

# Membuat model Multiple Linear Regression dari Training set
from sklearn.linear_model import LinearRegression
regressor = LinearRegression()
regressor.fit(X_train, y_train)

# Memprediksi hasil Test set
y_pred = regressor.predict(X_test)

# Memilih model multiple regresi yang paling baik dengan metode backward propagation
import statsmodels.api as sma
X = sma.add_constant(X)
import statsmodels.formula.api as sm
X_opt = X[:, [0, 1, 2, 3, 4, 5]]
regressor_OLS = sm.OLS(endog = y, exog = X_opt).fit()
regressor_OLS.summary()
X_opt = X[:, [0, 1, 3, 4, 5]]
regressor_OLS = sm.OLS(endog = y, exog = X_opt).fit()
regressor_OLS.summary()
X_opt = X[:, [0, 3, 4, 5]]
regressor_OLS = sm.OLS(endog = y, exog = X_opt).fit()
regressor_OLS.summary()
X_opt = X[:, [0, 3, 5]]
regressor_OLS = sm.OLS(endog = y, exog = X_opt).fit()
regressor_OLS.summary()
X_opt = X[:, [0, 3]]
regressor_OLS = sm.OLS(endog = y, exog = X_opt).fit()
regressor_OLS.summary()

Penjelasan:

  • Line 2 sampai line 4 mengimpor library yang diperlukan.
  • Line 7 mengimpor datasetnya.
  • Line 8 mendefinisikan variabel independennya (biaya R&D, marketing, administrasi dan wilayah). Di sini kita belum mendefinisikan dummy variabel.
  • Line 9 menampilkan X sebagai data frame. Jika Anda menggunakan spyder versi lama, maka ketika Anda double klik X di variable inspector, Anda bisa melihat isi datanya. Namun di spyder yang baru hal ini tidak dimungkinkan, karena kolom wilayah bukanlah tipe data numerik. Jika Anda ingin sekilas melihat data X, cukup ketik X di console lalu tekan ENTER. Namun Anda juga bisa mengeksekusi line 9 untuk melihat isinya sebagai sebuah data frame.
  • Line 10 mendefinisikan variabel dependennya.
  • Line 13 mengimpor class LabelEncoder dan OneHotEncoder untuk membuat dummy variabel.
  • Line 14 mempersiapkan objek labelencoder untuk proses dummy variabel berikutnya. Jika Anda lupa apa saja yang harus dituliskan, cukup copy paste dari data format di pembahasan preprocessing di link ini.
  • Line 15 menentukan kolom ke berapa yang ingin dibuat dummy variabelnya. Karena indeks python dimulai dari nol, maka kolom terakhir adalah kolom ketiga. Line 15 hanya mengonversi kategori menjadi 0, 1 atau 2.
  • Label 16 melakukan transformasi OneHotEncoder hasil dari line 15 menjadi dummy variabel yang kita inginkan (nilai 0 atau 1). Perlu diingat karena punya 3 jenis kategori, maka kita tuliskan categorical_features = [3].
  • Line 17 mengeksekusi proses transformasi kolom wilayah objek X menjadi dummy variabel. Anda sekarang bisa klik 2x pada variable inspector variabel X. Di situ bisa dilihat 3 kolom pertama sudah menjadi dummy variable. Namun perlu diingat kita hanya memerlukan n-1 dummy variable. Oleh karena itu line kita hilangkan salah satu kolom di line 18.
  • Line 18 menghilangkan salah satu kolom dummy variabel agar tidak terjebak dummy variable trap.

Hasil transformasi dataset menjadi 2 dummy variabel akan tampak sebagai berikut:

 

 

  • Line 23 mengimpor class dan sublibrary yang diperlukan untuk membagi dataset ke train set dan test set.
  • Line 24 membagi data 20% ke test set dan 80% ke train set.
  • Line 27 mengimpor class LinearRegression untuk membuat model regresi.
  • Line 28 mempersiapkan objek regressor sebagai model regresi kita nantinya.
  • Line 29 membuat model regresi dari dataset X_train dan y_train.
  • Line 32 memprediksi hasil test set, di mana kita ingin membandingkan nilai y_pred dengan y_test nantinya.

Perbandingan antara y_pred dengan y_test akan tampak sebagai berikut:

 

 

Bisa dilihat bahwa model kita cukup baik memprediksi nilai profit 50 perusahaan ini. Terlihat dari selisih yang tidak terlalu besar untuk setiap barisnya antara y_pred dan y_test. Walau demikian, tentunya tidak cukup jika menilai baik tidaknya hanya membandingkan secara visual saja. Kita memerlukan cara kuantitatif dengan menggunakan statistik.

Oleh karena itu, untuk menguji apakah model ini sudah baik, dan variabel independen mana saja yang signifikan atau tidak, kita melakukan metoda backpropagation atau backward propagation. Metode ini memulai dengan memasukkan semua variabel independen, kemudian melucuti satu per satu variabel independen yang tidak signifikan sampai ditemukan model yang terbaik (converged). Bagaimana memilih variabel independen yang dilucuti? Adalah variabel independen yang p value nya di atas 0.05. Kemudian kita uji lagi, jika ada variabel independen yang p value nya di atas 0.05 kita buang lagi. Begitu seterusnya sampai converged.

  • Line 35 mengimpor library stastmodels.api sebagai sma
  • Line 36 memerlukan penjelasan yang sedikit panjang. Berikut penjelasannya:

Seperti kita ketahui, fungsi umum dari multi regresi adalah: y = a_{{0}} + a_{{1}}*x_{{1}} + a_{{2}}*x_{{2}} + ... + a_{{n}}*x_{{n}}

Di situ kita memiliki konstanta a_{{0}}  di mana ia juga sebenarnya dikalikan dengan variabel x_{{0}}  dengan nilai 1. Library statsmodels tidak mempertimbangkan hal ini. Sementara library LinearRegression di line 27 sudah mempertimbangan variabel x_{{0}} . Oleh karena itu line 36 adalah memberikan tambahan kolom baru di objek X dengan nilai 1.

Jika dilihat objek X yang baru akan tampak sebagai berikut:

 

 

Kita lihat di kolom paling kiri merupakan nilai x_{{0}}  yang semua barisnya adalah 1.

  • Line 37 mengimpor library statsmodels.formula.api sebagai sm.
  • Line 38 membuat objek X_opt sebagai daftar semua variabel independen yang ingin dianalisis. Di sini kita memilih semua baris, dan semua kolom. Tentunya kita bisa menuliskannya dengan X[ : , : ]. Tapi kita menuliskannya dengan X[ : , [0,1,2,3,4,5]] untuk nantinya mempermudah kita membuang satu per satu variabel independen yang tidak signifikan.
  • Line 39 membuat objek X_OLS sebagai objek untuk melihat performa model kita dengan menggunakan metode OLS (Ordinary least squares) dari library sm kita tadi. Anda tentu bisa selalu menginspeksi parameter apa saja yang diperlukan dengan menekan CTRL+i di OLS(). Jika tidak muncul setelah klik CTRL+i, maka coba ketik sm.OLS? di console. Parameter yang diperlukan adalah endog sebagai variabel dependen, dan exog sebagai variabel independennya (dengan sudah menambahkan x_{{0}} ). Oleh karena itu endog = y, dan exog = X_opt. Jangan lupa menambahkan metode .fit() di belakangnya.

 

Tips: Dalam penulisan metode di python ada dua cara, langsung dan tidak langsung. Berikut penullisan yang berbeda namun hasilnya sama, tergantung bagaimana gaya dan kenyamanan Anda.

# STYLE 1
regressor_OLS = sm.OLS(endog = y, exog = X_opt).fit()
regressor_OLS.summary()

# STYLE 2
regressor_OLS = sm.OLS(endog = y, exog = X_opt)
hasil = regressor_OLS.fit()
hasil.summary()

 

  • Jika Anda eksekusi line 40, hasilnya akan tampak sebagai berikut:

 

 

Lihatlah pada bagian p value nya (P>|t|). Di mana semakin kecil nilai maka semakin baik variabel tersebut. Semakin tinggi semakin buruk. Cara mudahnya adalah nilai di atas 0.05 maka ia tidak signifikan sebagai sebuah variabel independen untuk memprediksi profit perusahaan.

Kita juga bisa melihat bahwa ada 5 variabel independen, dengan penjelasan sebagai berikut:

  • const adalah nilai koefisien untuk x_{{0}}
  • x1 dan x2 adalah dummy variabel untuk variabel independen Wilayah (p value x1 = 0.95, dan x2 = 0.99)
  • x3 adalah variabel independen R&D (p value x3 < 0.05)
  • x4 adalah variabel independen administrasi (p value x4 = 0.6)
  • x5 adalah variabel independen marketing (p value x5 = 0.12)

Setelah dilihat maka variabel independen dengan p value terbesar harus dibuang. Maka kita membuang x2.

Jika Anda eksekusi terus dari line 41 sampai line 52 maka variabel independen yang tersisa adalah variabel R&D dengan tampilan sebagai berikut:

 

 

Dengan demikian kita bisa menuliskan fungsi multiple regresinya sebagai berikut:

 

Profit = 49030 + 0.8543*BiayaR.D

 

Jika dirasa proses backpropagation ini terlalu melelahkan karena harus manual satu per satu mengeluarkan variabel independen yang tidak signifikan, maka Anda bisa melakukannya secara otomatis dengan beberapa lines yang sudah saya siapkan di bawah ini:

import statsmodels.formula.api as sm
def backwardElimination(x, sl):
    numVars = len(x[0])
    for i in range(0, numVars):
        regressor_OLS = sm.OLS(y, x).fit()
        maxVar = max(regressor_OLS.pvalues).astype(float)
        if maxVar &gt; sl:
            for j in range(0, numVars - i):
                if (regressor_OLS.pvalues[j].astype(float) == maxVar):
                    x = np.delete(x, j, 1)
    regressor_OLS.summary()
    return x

SL = 0.05
X_opt = X[:, [0, 1, 2, 3, 4, 5]]
X_Modeled = backwardElimination(X_opt, SL)

 

 

Bahasa R

# Mengimpor dataset
dataset = read.csv('50_Startups.csv')
 
# Encode categorical data
dataset$Wilayah = factor(dataset$Wilayah,
                       levels = c('New York', 'California', 'Florida'),
                       labels = c(1, 2, 3))
 
# Membagi data ke Training set and Test set
# install.packages('caTools')
library(caTools)
set.seed(123)
split = sample.split(dataset$Profit, SplitRatio = 0.8)
training_set = subset(dataset, split == TRUE)
test_set = subset(dataset, split == FALSE)
 
# Membangun model Multiple Linear Regression dari Training set
regressor = lm(Profit ~ R.D + Administrasi + Marketing + Wilayah, data = training_set)
summary(regressor)
 
# Alternatif penulisan
regressor2 = lm(Profit ~ ., data = training_set)
 
# Memprediksi hasil Test set
y_pred = predict(regressor, newdata = test_set)
 
# Membandingkan y_pred dengan y_test
perbandingan = data.frame(y_pred,test_set$Profit)
 
# Membangun model terbaik menggunakan metode Backpropagation atau Backward Elimination
regressor = lm(Profit ~ R.D + Administrasi + Marketing + Wilayah, data = dataset)
summary(regressor)
regressor = lm(Profit ~ R.D + Administrasi + Marketing + factor(Wilayah, exclude = 2), data = dataset)
summary(regressor)
regressor = lm(Profit ~ R.D + Administrasi + Marketing, data = dataset)
summary(regressor)
regressor = lm(Profit ~ R.D + Administrasi + Marketing, data = dataset)
summary(regressor)
regressor = lm(Profit ~ R.D, data = dataset)
summary(regressor)

Penjelasan:

  • Line 2 mengimpor dataset yang diperlukan.
  • Line 5 meng-encode data kategori Wilayah. Catatan: di R membuat dummy variabel tidak serumit di Python.
  • Line 11 mengimpor library untuk membagi data ke test dan training set.
  • Line 12 adalah menentukan angka bilangan random (bebas).
  • Line 13 membagi data menjadi dua dengan 80% untuk train set dan 20% test set dengan mengonversi ke tipe data boolean (TRUE or FALSE)
  • Line 14, jika bernilai TRUE, maka ia menjadi train set.
  • Line 15, jika bernilai FALSE, maka ia menjadi test set.
  • Line 18 membuat model regresi berdasarkan training set.
  • Line 19 menampilkan ringkasan statistik model regresi yang kita buat. Tampilannya akan tampak sebagai berikut:

 

 

Dapat dilhat, jika kita menggunakan R, maka R secara otomatis akan mempertimbangkan masalah jebakan dummy variabel, sehingga kita tidak perlu khawatir. R secara otomatis membuang 1 dari 3 kategori, maka yang ditampilkan dalam ringkasan statistiknya hanyalah wilayah2 dan wilayah3.

Pada kolom Pr(>|t|) dapat dilihat nilai p value dari masing-masing variabel. Semakin kecil (apalagi jika signifikan di bawah 0.05, ditandai dengan simbol asterik *), maka variabel tersebut secara signifikan berpengaruh terhadap variabel dependen. Sebaliknya, jika semakin besar maka semakin tidak signifikan (semakin kecil pengaruhnya terhadap variabel dependen).

Karena kita ingin memilih model regresi yang terbaik (paling bisa menjelaskan pengaruh antara variabel dependen dan independen), maka kita buang variabel dengan nilai p value terbesar yaitu wilayah2. Metode yang digunakan adalah backpropagation.

Pada summary(regressor) kita sebenarnya bisa menduga bahwa variabel R.D lah yang memiliki pengaruh paling kuat terhadap Profit perusahaan.

  • Line 22 adalah cara lain menuliskan fungsi regresi dalam R. Daripada menuliskan satu demi satu variabel independen, Anda bisa cukup menuliskan titik (.) setelah tanda ~.
  • Line 25 adalah memprediksi hasil y_pred melalui pembelajaran train_set. Nantinya kita akan bandingkan y_pred dengan y_test.
  • Line 28 membuat data frame gabungan antara y_pred dengan kolom profit dari test_set. Perbandingannya tampaksebagai berikut:

 

 

Secara visual, prediksi kita terhadap y_test cukup baik. Walau demikian, kita harus menentukan model terbaik kita secara statistik (kuantitatif).

  • Line 31 sekilas sama dengan line 18, namun berbeda.  Kali ini kita menggunakan dataset awal, karena ingin mencari model terbaik. Kita sudah tidak tertarik lagi membangun hubungan antara X_train dan y_train.
  • Line 32 memberikan ringkasan statistik dari line 31. Hasilnya tampak sebagai berikut:

 

 

Karena variabel wilayah2 memiliki p value terbesar, maka wilayah2 kita buang.

  • Line 33 kita membuang wilayah2 dari perbandingan model. Perhatikan bagaimana cara mengeluarkan wilayah2, yaitu dengan perintah factor(Wilayah, exclude = 2)
  • Line 35 sampai line 40 didapat bahwa variabel independen yang tersisa adalah R.D, sehingga hasilnya nampak sebagai berikut:

 

 

Jika kita tulis fungsi regresinya maka tampak sebagai berikut:

 

Profit = 49030 + 0.8543*BiayaR.D

 

Fungsinya sama persis antara python dan R, sehingga Anda bisa memilih mana di antara keduanya yang paling nyaman.

Jika Anda merasa metode backward propagation terlalu melelahkan, karena harus satu demi satu kita mengeluarkan variabel independen yang tidak signifikan, maka saya siapkan script R di bawah ini yang secara otomatis menyaring variabel independen yang tidak signifikan tersebut. Sehingga hidup Anda terasa lebih mudah.

backwardElimination &lt;- function(x, sl) {
    numVars = length(x)
    for (i in c(1:numVars)){
        regressor = lm(formula = Profit ~ ., data = x)
        maxVar = max(coef(summary(regressor))[c(2:numVars), "Pr(&gt;|t|)"])
        if (maxVar &gt; sl){
            j = which(coef(summary(regressor))[c(2:numVars), "Pr(&gt;|t|)"] == maxVar)
            x = x[, -j]
        }
        numVars = numVars - 1
    }
    return(summary(regressor))
}

SL = 0.05
dataset = dataset[, c(1,2,3,4,5)]
backwardElimination(training_set, SL)

Perlu diingat! Kebetulan model multiple regresi kita hanya menyisakan 1 variabel independen (yaitu R&D), ini hanyalah kebetulan saja. Bisa jadi untuk kasus lain, semua variabel independen menjadi signifikan, atau bisa jadi hanya tersisa dua, tiga dan seterusnya.

Kesimpulan apa yang bisa diambil?

Sampai di sini kita bisa mengerti bahwa profit 50 perusahaan sangat ditentukan oleh biaya R&D nya, sehingga semakin banyak sebuah perusahaan berinvestasi di biaya pengembangan ini, maka semakin tinggi profitnya. Selian itu kita juga mengetahui bahwa biaya marketing, administrasi dan jenis wilayah tidak berdampak signifikan terhadap profit perusahaan.

Apakah hasil ini valid (sesuai) bagi semua perusahaan (perusahaan lain di luar 50 perusahaan ini)? Tentu saja tidak, karena model ini dibangun berdasarkan data 50 perusahaan yang kita miliki. Tapi setidaknya hasil dari model regresi ini bisa menjadi senjata untuk membuat strategi tentang bagaimana cara meningkatkan profit ke depannya.

Saya harap Anda bisa membuat model multiple regresi machine learning Anda sendiri. Anda bisa mengaplikasikannya untuk berbagai persoalan real. Misal, apakah ada hubungan antara nilai ujian dengan usia pelajar, jenis kelamin, lama belajar, uang bulanan dari orang tua, uang spp yang dibayarkan, dan faktor-faktor lainnya.

Multiple regresi juga bisa diaplikasikan untuk data biologis seperti data pertumbuhan tanaman, data DNA, data medis, data eskerimen dan masih banyak lagi. Anda bisa memasukkan banyak variabel independen yang Anda duga berpengaruh terhadap variabel dependen Anda.

Jika ada pertanyaan, kita bisa berdiskusi di bagian komentar.

Semoga bermanfaat.

Leave a Reply

avatar
  Subscribe  
Notify of