Email : je070601@gmail.com
Instagram : https://www.instagram.com/marvis.zerex/
RPubs : https://rpubs.com/invokerarts/
Linkedin : https://www.linkedin.com/in/jeffry-wijaya-087a191b5/
Majors : Business Statistics
Address : ARA Center, Matana University Tower
Jl. CBD Barat Kav, RT.1, Curug Sangereng, Kelapa Dua, Tangerang, Banten 15810.
Seperti yang Saya sudah katakan sebelumnya, Model Regresi digunakan untuk menunjukan hubungan antar variabel dengan garis yang memnuhi data observasi. Model regresi linear menggunakan garis lurus, sedangkan model regresi logistik dan linear menggunakan garis melengkung. Regresi memungkinkan Anda untuk memperkirakan variabel tetap berubah seperti variabel bebasnya berubah.
Dalam bagian ini, kita akan fokus dalam bagaimana menggunakan Regresi Linear Sederhana untuk memperkirakan hubungan antara 2 variabel-variabel kuantitatif. Kita juga dapat mengugnakan regresi linear sederhana untuk mengetahui:
Contoh:
Asumsikan Anda adalah peneliti sosial yang tertarik dengan hubungan antara pendapatan dan kesenangan pelanggan. Anda melakukan survei kepada 500 orang yang berpenghasilan dengan interval dari $15k hingga $75k dan bertanya kepada mereka tingkat kebahagiaan dari skala 1 sampai 10. Kedua variabel Anda baik variabel tetap (pendapatan) dan variabel bebas (kesenangan) bersifat kuantitatif, jadi anda dapat melakukan analisis regresi untuk melihat apakah ada hubungan linear diantara kedua variabel tersebut.
Sebuah model regresi linear menasumsikan bahwa ada hubungan linear diantara ekspetasi kondisi dari kedua variabel kuantitatif:
where,
Data untuk contoh ini terdapat pada package R yaitu PoEdata. Pertama-tama, Anda harus menginstall package “PoEdata”, ketik perintah berikut pada bagian Console Rstudio Anda:
install.packages("devtools") # Jika belum di install
library(devtools) # Meng-load packagenya
install_git("https://github.com/ccolonescu/PoEdata") # Mengintall `PoEdata`Kemudian, pertimbangkan data set:
library(PoEdata) # Meng-load package `PoEdata`
library(DT) # Menggunakan fungsi datatable
data("food", package="PoEdata")
datatable(food)Akan selalu menjadi ide yang baik untuk memeriksa data secara visual menggunakan diagram scatter, dimana dapat dibuat menggunakan fungsi plot(). Gambar dibawah adalah diagram scatter dari pengeluaran makanan berdasarkan penghasilan, dapat dilihat bahwa ada hubungan positif antara pendapatan dan pengeluaran makanan.
plot(food$income, food$food_exp,
ylim=c(0, max(food$food_exp)),
xlim=c(0, max(food$income)),
xlab="weekly income in $100",
ylab="weekly food expenditure in $",
type = "p")Untuk data pengeluaran makanan, model regresinya akan menjadi
\[ \begin{align} \tag{1} Y&=\beta_1+\beta_2 *X+e \\ food_{exp}&=\beta_1+\beta_2*income+e \end{align} \]
Fungsi R untuk memperkirakan model regresi linear adalah lm(y~x, data) dimana, jika digunakan dengan sendirinya tidak menunjukan hasil yang jelas; Sangat berguna untuk memberikan nama tersebut model sepertimod1, kemudian menunjukan hasilnya menggunakan summary(mod1).
mod1 <- lm(food_exp ~ income, data = food) # Model Linear
smod1<- summary(mod1) # Mendapatkan Summary
smod1##
## Call:
## lm(formula = food_exp ~ income, data = food)
##
## Residuals:
## Min 1Q Median 3Q Max
## -223.025 -50.816 -6.324 67.879 212.044
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 83.416 43.410 1.922 0.0622 .
## income 10.210 2.093 4.877 1.95e-05 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 89.52 on 38 degrees of freedom
## Multiple R-squared: 0.385, Adjusted R-squared: 0.3688
## F-statistic: 23.79 on 1 and 38 DF, p-value: 1.946e-05
Jika anda tertarik hanya dengan beberapa hasil dari regresi, seperti estimasi koefisien, Anda dapat mendapatkannya dengan menggunakan perintah khusus. seperti coef().
## [1] 83.416
## [1] 10.20964
Daru hasil tersebut, kita dapat menyatakan persamaan estimasi regresi linear kita menjadi:
\[ \begin{align} \tag{2} \hat{Y}&=b_1+b_2*X \\ \hat{food_{exp}}&=83.416+10.20964*income \end{align} \]
Tingkat intercept yaitu, \(\beta_1\) , biasanya hanya sedikit penting dalam model ekonometrik; kita akan sering lebih tertarik pada tingkat slope yaitu, \(\beta_2\). Nilai perkitaan dari \(\beta_2\) menunjukan bahwa pengeluaran makanan untuk kebanyakan keluarga meningkan sebesar 10.209643 ketika penghasilan keluarga naik sebesar 1, dimana dalam kasus ini adalah $100. Fungis R abline() menambahkan garis regresi pada diagram scatter yang sudah dibuat sebelumnya menjadi seperti dibawah ini:
plot(food$income, food$food_exp,
xlab="weekly income in $100",
ylab="weekly food expenditure in $",
type = "p")
abline(b1,b2)Bagaimana cara mengambil berbagai kesimpulan regresi? Untuk mendapatkan sebuah kesimpukan, Anda cukup mengetik nama model Anda, diikuti dengan tanda $ dan nama dari hasil yang ingin anda cari. Sebagai contoh, juga kita ingin mengetahui vektor koefisien dari mod1, kita menulisnya sebagaimod1\$coefficients dan smod1\$coefficients:
## (Intercept) income
## 83.41600 10.20964
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 83.41600 43.410163 1.921578 6.218242e-02
## income 10.20964 2.093264 4.877381 1.945862e-05
Seperti yang kita lihat sebelumnya, beberapa hasil dapat dicari menggunakan fungsi khusus seperti coef(mod1), resid(mod1), fitted(mod1), and vcov(mod1).
Perkiraan tingkat regresi, \(b_1\) dan \(b_2\) memungkinkan kita untuk memprediksi ekspetasi pengeluaran makanan untuk setiap penghasilan yang diketahui. Semua yang kita butuhkan hanya cukup memasukan estimasi nilai parameter dan penghasilan yang ada pada persamaan seperti persamaan ke-2. Sebagai contoh, nilai ekspetasi dari food_exp Untuk penghasilan $2000 dihitunga dalam persamaan ke-3 . (Jangan lupa untuk membagi nilai penghasilan dengan 100 karena data variabel yang kita gunakan dihitung dalam ratusan dolar)
\[ \begin{align} \tag{3} \hat{food_{exp}}&=83.416+10.20964*20\\ & = $ 287.608861 \end{align} \]
Namun, R akan melakukan kalkulasi dengan fungsi predict(). Mari kita coba memperluas contoh yang kita punya dengan menggunakan lebih dari satu penghasilan untuk memrediksi pengeluaran makanan, anggap penghasilan = $2000, $2500, dan $2700. Fungsi predict() dalam R membutuhkan nilai baru dari variabel bebas dibentuk dalam suatu data frame. Bahkan ketika kita hanya ingin memprediksi 1 penghasilan, kita membutuhkan struktur data framte yang sama. Dalam R, sebuah kelompok angka dinyanakan dalam struktur c(). Perintah berikut dapat mempermudah anda mengerti.
mod1 <- lm(food_exp~income, data=food)
newx <- data.frame(income = c(20, 25, 27))
yhat <- predict(mod1, newx)
yhat # menunjukan hasil## 1 2 3
## 287.6089 338.6571 359.0764
Koefisien regresi \(b_1\) dan \(b_2\) adalah variabel random, karena mereka bergantung kepada sampel. Dalam kasus ini, kita mengulang sampel pada nilai koefisien regresi. Kita akan membuat angka dari subsample random dari data food dan menghitung ulang \(b_1\) and \(b_2\). Subsample random dapat dibuat menggunakan fungsi sample(), contoh begitu mengilustrasikan hanya untuk mencari \(b_2\).
N <- nrow(food) # mengetahui jumlah data observasi dalam dataset
C <- 50 # jumlah subsample yang diharapkan
S <- 38 # jumlah sample yang diharapkan
sumb2 <- 0
for (i in 1:C){ # pengulangan jumlah dari subsample
set.seed(3*i) # kode yang berbeda untuk setiap subsample
subsample <- food[sample(1:N, size=S, replace=TRUE), ]
mod2 <- lm(food_exp~income, data=subsample)
#jumlah b2 untuk semua subsample:
sumb2 <- sumb2 + coef(mod2)[[2]]
}
print(sumb2/C, digits = 3)## [1] 9.89
Hasilnya, \(b_2=9.89\), adalah rata-rata dari estimasi 50 \(b_2\).
Banyak proses yang diperlukan untuk mengestimasi varians dan kovarians dari sebuah koefisien regresi. R menyatukannya dengan fungsi vcov():
## [1] 1884.442
## [1] 4.381752
## [1] -85.90316
Terkadang, diagram scatter atau beberapa perkiraan secara teori menunjukan hubungan non-linear. Hubungan non-linear yang paling sering ditemukan ada 2, yaitu yang berhubungan dengan logaritma dari variabel tetap atau bebas dan persamaan polinomial. Model kuadrat membutuhkan kuadrat dari variabel bebas.
\[
\begin{align}
\tag{4}
Y=\beta_1+\beta_2*X^2+e\\
\end{align}
\] Dalam R, variabel bebas melibatkan operasi matematika yang meliputi fungsi regresi dengan fungsi I(). Contoh berukut menggunakan dataset br dari packkage PoEdata, yang mana mencakup harga penjualan dan atribut dari 1080 rumah di Baton Rouge, LA. Harga dalam penjualan dihitung dalam dolar, dan dqft merupakan luas tanah dalam satuan kaki kuadrat.
data(br)
mod2 <- lm(price~I(sqft^2), data=br)
b1 <- coef(mod2)[[1]]
b2 <- coef(mod2)[[2]]
sqftx=c(2000, 4000, 6000) # nilai yang diberikan untuk luas tanah
pricex=b1+b2*sqftx^2 # harga yang sesuai dengan luas tanah yang diberikan
DpriceDsqft <- 2*b2*sqftx # efek marginal dalam luas tanah pada harga
elasticity=DpriceDsqft*sqftx/pricex
b1; b2; DpriceDsqft; elasticity # menunjukan hasil## [1] 55776.57
## [1] 0.0154213
## [1] 61.68521 123.37041 185.05562
## [1] 1.050303 1.631251 1.817408
Sekarang kita ingin menunjukan diagram scatter dan melihat bagaimana persamaan kuadrat mencocokan data. R Chuck selanjutnya menunjukan dua cara untuk membuat grafik tersebut. Pertama-tama, cukup menunjukan persamaan kuadrat pada diagram scatter, menggunakan fungsi curve(); kedua, gunakan garis persamaan, yang mana dibutuhkan pengurtan data secata menaik berdasarkan sqft sebelum model regresi dievaluasi
plot(br$sqft, br$price, xlab="Total square feet",
ylab="Sale price, $", col="grey")
#menambahkan persamaan kuadrat pada scatterplot:
curve(b1+b2*x^2, col="red", add=TRUE) Cara alternatif untuk menggambar kurva yang cocok adalah:
ordat <- br[order(br$sqft), ] #mengurutkan data berdasarkan `sqft`
mod2 <- lm(price~I(sqft^2), data=ordat)
plot(br$sqft, br$price,
main="Dataset ordered after 'sqft' ",
xlab="Total square feet",
ylab="Sale price, $", col="grey")
lines(fitted(mod31)~ordat$sqft, col="red")Model log-linear meregresi log dari variabel tetap pada hasil linear dari variabel bebas. (kecuali sudah ditentukan, notasi log menunjuk kepada logaritma natural, mengikuti konvensi umum pada bidang ekonomi):
\[ \begin{align} \tag{5} Log(Y)=\beta_1+\beta_2*X+e\\ \end{align} \]
Salah satu alasan menggunakan logaritma pada variabel bebas adalah membuat distribusinya mendekat kepada distribusi normal. Mari kita tunjukan histogram dari harga dan log(harga) untuk membandingkannya
Perbedaan setelah dilogaritmakan dengan sebelum dilogaritmakan adalah * Distribusinya, dengan menggunakan logaritma, kita mendekatkan distribusinya kepada distribusi normal. * Outliers-nya, dengan menggunakan logaritma, kita mengurangi banyak outliers yang ada pada kumpulan data tersebut.
Kita tertarik dalam mengestimasi koefisien dan intercept mereka, dalam nilai harga yang sesuai dan efek marginal dari peningkatan luas tanah berdasarkan harga.
## (Intercept) sqft
## 1.083860e+01 4.112689e-04
Koefisiennya adalah \(b_1=10.84\) dan \(b_2=0.00041\), Menunjukan bahwa peningkatan luas tanah (sqft) dari apartemen setiap 1 unit (1 sqft) meningkatkan harga apartementnta sebesar 0.041 persen. Jadi, untuk rumah yang berharga $100,000, setiap peningkatan 100 sqft akan meningkatkan harga sekitar 100∗0.041 persen, yang bernilai $4112.7. Secara general, efek marginal dari peningkatan \(X\) pada \(Y\) dalam persamaan 5 adalah
\[ \begin{align} \tag{6} {dy\over dx}=\beta_2 * Y \end{align} \]
Dan untuk elastisitasnya adalah
\[ \begin{align} \tag{7} e= {dy\over dx} {X\over Y}=\beta_2 * X \end{align} \]
Beberapa perintah dibawah, menunjukan bagaimana cara menggambar nilai kurva dari model log-linear dan bagaimana cara menghitung efek marginal dan elastisitas untuk nilai tengah harga pada dataset. Nilai yang tepat dapat dihitung menggunakan formula berikut:
\[ \begin{align} \tag{8} \hat{Y}= e^{b_1+b_2*X} \end{align} \]
ordat <- br[order(br$sqft), ] #menurutkan dataset
mod4 <- lm(log(price)~sqft, data=ordat)
plot(br$sqft, br$price, col="grey")
lines(exp(fitted(mod4))~ordat$sqft,
col="blue", main="Log-linear Model")pricex<- median(br$price)
sqftx <- (log(pricex)-coef(mod4)[[1]])/coef(mod4)[[2]]
(DyDx <- pricex*coef(mod4)[[2]])## [1] 53.46495
## [1] 0.9366934
R Memungkinkan kuta untuk menghitung kuantitas yang sama untuk beberapa (luas tanah, harga) pasangan dalam waktu yang sama, sepetri yang dilihatkan berikut:
b1 <- coef(mod4)[[1]]
b2 <- coef(mod4)[[2]]
#mengambul beberapa nilai untuk luas tanah(sqft):
sqftx <- c(2000, 3000, 4000)
#mengestikasi harga untuk luas tanah dan menambahkan lagi:
pricex <- c(100000, exp(b1+b2*sqftx))
#menghitung ulang luas tanah (sqft) untuk seluruh harga:
sqftx <- (log(pricex)-b1)/b2
#menghitung kan menununjukan elastisitasnya:
(elasticities <- b2*sqftx) ## [1] 0.6743291 0.8225377 1.2338066 1.6450754
Indikator atau variabel biner menunjukan wujud dari beberapa atribut dari unit yang di observasi, seperti jenis kelamin maupun suku jika unit yang di observasi bersifat individual atau lokasi jika unit yang di observasi adalah rumah. Dalam dataset utown, variabel utown adalah 1 jika rumah dekat dengan kampus dan 0 untuk sebaliknya. Ini adalah contoh model regresi linear yang melibatkan variabel utown:
\[ \begin{align} \tag{9} price_i= \beta_1+\beta_2*utwon_i \end{align} \]
Koefisien dari variabel dalam model linear sederhana sama dengan perbedaan rata-rata harga dari kedua kategori; koefisien intercept dari model dalam persamaan 9 sama dengan rata-rata harga dari rumah yang tidak dekat dengan kampus. Mari pertama-tama kita hitung rata-rata dari harga untuk setiap kategori, dimana dinyatakan dengan kode price0bar dan price1bar:
data(utown)
price0bar <- mean(utown$price[which(utown$utown==0)])
price1bar <- mean(utown$price[which(utown$utown==1)])Hasil yang didapatkan adalah: \(\bar{price}=\) 277.2416012 dekat dengan kampus, dan \(\bar{price}=\) 215.7324948 untuk yang tidak dekat. Sekarang kita menunjukan bahwa hasil yang sama dihasilkan dari koefisien model regresi dalam Persamaan 9 :
## [1] 215.7325
## [1] 61.50911
Hasil yang didapatkan adalah: \(\bar{price}=b_1=\) 215.7324948 untuk yang rumahnya jauh dari kampus, dan \(\bar{price}=b_1+b_2=\) 277.2416012 yang rumahnya dekat dari kampus.
Simulasi Monte Carlo membuat nilai acak untuk variabel tetap ketika koefisien regresi dan distribusi random diberikan. Contoh berikut dilakukan untuk menentukan distribusi dari variabel bebas pada model pengeluaran makanan dalam persamaan 2,
N <- 40
x1 <- 10
x2 <- 20
b1 <- 100
b2 <- 10
mu <- 0
sig2e <- 2500
sde <- sqrt(sig2e)
yhat1 <- b1+b2*x1
yhat2 <- b1+b2*x2
curve(dnorm(x, mean=yhat1, sd=sde), 0, 500, col="blue")
curve(dnorm(x, yhat2, sde), 0,500, add=TRUE, col="red")
abline(v=yhat1, col="blue", lty=2)
abline(v=yhat2, col="red", lty=2)
legend("topright", legend=c("f(y|x=10)",
"f(y|x=20)"), lty=1,
col=c("blue", "red"))Kemudian, kita menghitung varians dari \(b_2\) dan menunjukan plot dari density funtion berikut.
\[ \begin{align} \tag{10} var(b_2)={\sigma^2 \over \sum(x_i−\bar{x})} \end{align} \]
x <- c(rep(x1, N/2), rep(x2,N/2))
xbar <- mean(x)
sumx2 <- sum((x-xbar)^2)
varb2 <- sig2e/sumx2
sdb2 <- sqrt(varb2)
leftlim <- b2-3*sdb2
rightlim <- b2+3*sdb2
curve(dnorm(x, mean=b2, sd=sdb2), leftlim, rightlim)
abline(v=b2, lty=2)Sekarang, dengan nilai \(b_1\), \(b_2\), dan error standar deviasi yang sama, kita dapat membuat set nilai untuk \(Y\) , regresi \(Y\) terhadap \(X\), dan menghitung nilai estimasi untuk koefisien \(b_2\)dan standar errornya.
set.seed(12345)
y <- b1+b2*x+rnorm(N, mean=0, sd=sde)
mod6 <- lm(y~x)
b1hat <- coef(mod6)[[1]]
b2hat <- coef(mod6)[[2]]
mod6summary <- summary(mod6) #summarynya terdapat nilai standard errors
seb2hat <- coef(mod6summary)[2,2]Hasilnya adalah \(b_2=11.64\) dan \(s_e(b_2)=1.64\). Kelebihan dari Simulasi Monte Carlo adalah, kemungkinan dari pengulangan estimasi tingkat regresi untuk angka yang besar akan secara otomatis membuat sample. Jadi, kida dapat mendapatkan angka yang besar untuk parameter, katakanlah \(b_2\), dan kemudian memperkitakan karakteristik samplingnya. Sebagai contoh, jika rata-rata dari nilai tersebut mendekati nilai asumsi yaitu \(b_2= 10\), kita dapat menyimpulkan bahwa estimator (metode dari mengestimasi parameter) adalah tidak bias.
Sekarang kita akan menggunakan nilai dari \(X\) dalam dataset food, dan menghasilkan \(Y\) menggunakan model linear dengan \(b_1=100\) dan \(b_2= 10.\)
data("food")
N <- 40
sde <- 50
x <- food$income
nrsim <- 1000
b1 <- 100
b2 <- 10
vb2 <- numeric(nrsim) #menyimpan estimasi nilai b2
for (i in 1:nrsim){
set.seed(12345+10*i)
y <- b1+b2*x+rnorm(N, mean=0, sd=sde)
mod7 <- lm(y~x)
vb2[i] <- coef(mod7)[[2]]
}
mb2 <- mean(vb2)
seb2 <- sd(vb2)Rata-rata dan standar deviasi dari estimasi 40 nilai dari \(b_2\) adalah, berturut-turut, 9.974985 dan 1.152632. Gambar dibawah menunjukan distribusi data yang telah disimulasikan dari \(b_2\) dengan yang dijelaskan secara teori.
plot(density(vb2))
curve(dnorm(x, mb2, seb2), col="red", add=TRUE)
legend("topright", legend=c("true", "simulated"),
lty=1, col=c("red", "black"))