Ekonometrik
Regresi Linier Sederhana
Putri Angelina Windjaya - 20194920010
February 09, 2021
Email : putriangelina865@gmail.com
Instagram : https://www.instagram.com/putriangelinaw
RPubs : https://rpubs.com/putriangelinaw/
1 Pengenalan
Seperti yang sudah saya sebutkan sebelumnya, Model Regresi digunakan untuk menjelaskan hubungan antara variabel dengan menyesuaikan garis dengan data observasi. MOdel regersi linier menggunakan garis lurus, sementara model regresi logistic dan nonlinier menggunakan garis kurva. Regresi dapat mengestimasi bagaimana variabel dependen berubah seiring variabel independen berubah.
Pada bagian ini kita akan fokus bagaimana menggunakan Regresi Linier Sederhana untuk mengestimasi hubungan antara dua variabel. Kita bisa juga menggunakan regresi linier sederhana untuk mengetahui:
- Seberapa kuat hubugnan antara dua variabel.
- Nilai dari variabel dependen pada nilai variabel independen tertentu.
Contoh:
Misalkan kamu adalah peneliti sosial yang tertarik dengan hubungan antara pendapatan dan kesenangan. Kamu survei 500 orang yang range pendapatannya dari $15k hingga $75k dan bertanya tentang kesenangan mereka dari 1 sampai 10. Variabel independen (pendapatan) dan variabel dependen (kesenangan) keduanya kuantitatif, jadi kamu bisa melakukan analisis regresi untuk melihat apakah terdapat hubungan linier antara keduanya.
2 Model Umum
Model regresi linier sederhana berasumsi bahwa terdapat hubungan linier antara ekspektasi bersyarat dari dua variabel kuantitatif:
dimana,
- \(Y\) sebagai prediktor, penjelasan, atau variabel dependen.
- \(X\) sebagai respoe, hasil, or variabel independen.
- \(\beta_1\) adalah parameter intercept atau koefisien
- \(\beta_2\) adalah parameter slope atau koefisien
- \(\hat{Y}\) adalah nilai estimasi \(Y\) dari \(X\)
- \(b_1\) dan \(b_2\) adalah nilai estimasi koefisien
3 Regresi Linier Sederhana pada R
Data yang akan digunakan adalah data dari R package PoEdata. Pertama, kamu harus install package PoEdata, dengan mengikuti koding dibawah ini pada console RStudio:
install.packages("devtools") # jika belum di install
library(devtools) # jalankan library
install_git("https://github.com/ccolonescu/PoEdata") # install `PoEdata`kemudian, meninjau dataset:
Dengan memvisualisasikan data dalam diagram scatter selalu menjadi ide yang bagus, yang mana bisa dibuat dengan function plot(). Gambar dibawah ini adalah diagram scatter dari biaya bahan pokok pada pendapatan, asusmikan bahwa terdapat hubungan positif antara pendapatan dengan biaya bahan 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")
3.1 Estimasi Regresi Linier
Untuk data biaya bahan makanan, model regresinya adalah sebagai berikut
\[ \begin{align} \tag{1} Y&=\beta_1+\beta_2 *X+e \\ food_{exp}&=\beta_1+\beta_2*income+e \end{align} \]
Function R untuk mengestimasi model regresi linier adalah lm(y~x, data) yang mana, digunakan dengan sendirinya tidak menunjukkan output apa pun; Itu berguna untuk memberikan nama pada model, misal mod1, lalu tampilkan hasil menggunakan summary(mod1).
mod1 <- lm(food_exp ~ income, data = food) # model linier
smod1<- summary(mod1) # hasil 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-05Jika kamu tertarik hanya pada beberapa hasil regresi, seperti koefisien yang diperkirakan, kamu bisa mengambilnya menggunakan function tertentu, seperti function coef().
## [1] 83.416## [1] 10.20964Dari hasil tersebut, kita bisa tulis estimasi regresi linier dengan:
\[ \begin{align} \tag{2} \hat{Y}&=b_1+b_2*X \\ \hat{food_{exp}}&=83.416+10.20964*income \end{align} \]
Parameter intercept, \(\beta_1\) , biasanya tidak terlalu penting dalam model ekonometrik; kita lebih tertarik dengan parameter slope, \(\beta_2\). Nilai estimasi \(\beta_2\) menunjukkan bahwa pengeluaran makanan untuk rata-rata keluarga meningkat sebesar 10.209643 ketika pendapatan keluarga meningkat 1 unit, yang dalam hal ini adalah $ 100. R function abline() menambahkan garis regression ke diagram scatter yang sebelumnya diplot:
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 hasil regresi? Untuk mengambil hasil tertentu kamu cukup merujuknya dengan nama objek, diikuti dengan tanda $ dan nama hasil yang ingin kamu ambil. Misalnya, jika kita menginginkan vektor koefisien dari mod1, kita menyebutnya 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-05Seperti yang sudah dilihat sebelumnya, akan tetapi, beberapa hasil ini bisa didapat menggunakan function yang lebih spesifik, seperti coef(mod1), resid(mod1), fitted(mod1), and vcov(mod1).
3.2 Prediksi Dengan Model Regresi Linier Sederhana
Parameter regresi yang diestimasi, \(b_1\) dan \(b_2\) memungkinkan kita untuk memprediksi biaya bahan makanan yang diharapkan untuk pendapatan tertentu. Yang perlu kita lakukan adalah memasukkan nilai parameter estimasi dan pendapatan yang diberikan ke dalam persamaan seperti Persamaan 2. Misalnya, nilai yang diharapkan dari food_exp untuk pendapatan $ 2000 dihitung dalam Persamaan 3. (Ingatlah untuk membagi pendapatan dengan 100, karena data untuk pendapatan variabel dalam ratusan dolar.)
\[ \begin{align} \tag{3} \hat{food_{exp}}&=83.416+10.20964*20\\ & = $ 287.608861 \end{align} \]
Bagaimanapun R melakukan perhitungan ini untuk kita dengan function predict(). Ayo kita sedikit memperluas contoh ke lebih dari satu pendapatan yang kita prediksi biaya bahan makanannya, katakanlah pendapatan = $ 2000, $ 2500, dan $ 2700. Function predict() di R mengharuskan nilai baru dari variabel independen diatur dalam bentuk tertentu, yang disebut data frame. Bahkan ketika kita hanya ingin memprediksi untuk satu pendapatan, kita membutuhkan struktur data frame yang sama. Di R, satu set bilangan disatukan menggunakan struktur c().
mod1 <- lm(food_exp~income, data=food)
newx <- data.frame(income = c(20, 25, 27))
yhat <- predict(mod1, newx)
yhat # hasilnya dicetak## 1 2 3
## 287.6089 338.6571 359.07643.3 Nilai Koefisien Regresi
Koefisien regresi \(b_1\) dan \(b_2\) adalah variabel random, karena mereka tergantung terhadap sampel. Pada kasus ini, kita mengulang sampel untuk menilai koefisien regresi. Ayo kita mengambil angka subsampel secara random dari data food dan menghitung ulang \(b_1\) dan \(b_2\) . subsampel yang random dapat dihasilkan menggunakan function sample(), contoh dibawah ini menampilkan hanya untuk \(b_2\).
N <- nrow(food) # banyak observasi dalam dataset
C <- 50 # memilih angka subsampel
S <- 38 # memilih ukuran sampel
sumb2 <- 0
for (i in 1:C){ # pengulangan di angka subsampel
set.seed(3*i)
subsample <- food[sample(1:N, size=S, replace=TRUE),]
mod2 <- lm(food_exp~income, data=subsample)
sumb2 <- sumb2 + coef(mod2)[[2]]}
print(sumb2/C, digits = 3)## [1] 9.89Hasilnya, \(b_2=9.88\), adalah rata-rata dari 50 estimasi \(b_2\).
3.4 Estimasi Var dan Cov
Banyak aplikasi memerlukan estimasi varians dan covarians dari koefisien regresi. R menyimpannya pada matriks vcov():
## [1] 1884.442## [1] 4.381752## [1] -85.903163.5 Arti dari Var dan Cov
Perhatikan gambar berikut:
Dari ilustrasi diatas kita dapat lebih mudah memahami arti varians ketika tinggi dan rendah. Misalkan ilustrasi diatas adalah sebuah papan pahanan atau biasa dikenal dart (dalam bahasa inggris) dan lingkaran yang hitam ditengah merupakan target kita. Ingat! asal kata varians adalah variasi, jadi ketika variasinya tinggi maka nilai-nilai data observasi memiliki perbedaan yang tinggi atau relatif berbeda, sedangkan jika variasinya rendah maka nilai-nilai data observasi relatif sama. Lalu coba kita perhatikan ketika bias nya rendah dan tinggi. Ketika bias nya tinggi maka data observasi berada jauh dari ekspektasi atau target kita, sedangkan ketika bias nya rendah maka data observasi berada dekat dengan ekspektasi atau target kita. Varians yang tinggi dalam varb1 mengartikan bahwa biaya bahan makanan dari pendapatan tiap sampel ini relatif berbeda. Dan varb2 itu rendah artinya seiring pendapatannya meningkat biaya bahan makanan akan relatif sama.
Nilai covb1b2 itu negatif mengartikan bahwa nilai pendapatan yang besar berkoresponden dengan nilai biaya bahan makanan yang kecil.
4 Hubungan Non-Linier
Terkadang scatter plot atau beberapa pertimbangan teoritis menunjukkan hubungan non-linier. Hubungan non-linier yang paling populer melibatkan logaritma dari variabel dependen atau independen dan fungsi polinomial. Model kuadrat membutuhkan kuadrat dari variabel independen.
\[
\begin{align}
\tag{4}
Y=\beta_1+\beta_2*X^2+e\\
\end{align}
\] Dalam R, variabel independen yang melibatkan operasi matematika dapat dimasukkan dalam persamaan regresi dengan function I(). Contoh dibawah ini menggunakan dataset br dari package PoEdata, yang mana terdapat harga jual dan lainnya dari 1080 rumah dalam Baton Rouge, LA. price adalah harga jual dalam dolar, dan sqft adalah luas permukaan dalam square feet.
data("br", package="PoEdata")
mod2 <- lm(price~I(sqft), data=br)
b1 <- coef(mod2)[[1]]
b2 <- coef(mod2)[[2]]
sqftx=c(2000, 4000, 6000) # beri nilai untuk sqft
pricex=b1+b2*sqftx^2 # harga untuk setiap sqft
DpriceDsqft <- 2*b2*sqftx # efek marginal dari sqft pada price
elasticity=DpriceDsqft*sqftx/pricex
b1; b2; DpriceDsqft; elasticity # cetak hasil## [1] -60861.46## [1] 92.74737## [1] 370989.5 741979.0 1112968.5## [1] 2.000328 2.000082 2.000036Sekarang kita ingin menggambar diagram scatter dan melihat bagaimana fungsi kuadrat cocok dengan data. Chunk berikutnya memberikan dua alternatif untuk membuat grafik. Yang pertama hanya menggambar fungsi kuadrat pada diagram scatter, menggunakan R function curve (); yang kedua menggunakan garis fungsi, yang memerlukan pengurutan kumpulan data dalam peningkatan nilai sqft sebelum model regresi dievaluasi, sehingga nilai pas yang dihasilkan juga akan keluar dalam urutan yang sama.
plot(br$sqft, br$price, xlab="Total square feet",
ylab="Sale price, $", col="grey")
#menambah kurva kuadrat kedalam scatter plot:
curve(b1+b2*x^2, col="red", add=TRUE) 
Cara alternatif untuk menggambar fitted curve:
ordat <- br[order(br$sqft),] #urutkan dataset setelah `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 dependen pada ekspresi linier dari variabel independen (kecuali ditentukan, notasi log adalah singkatan dari logaritma natural, mengikuti konvensi umum di bidang ekonomi):
\[ \begin{align} \tag{5} Log(Y)=\beta_1+\beta_2*X+e\\ \end{align} \]
Salah satu alasan untuk menggunakan log dari variabel independen adalah agar distribusinya mendekati distribusi normal. Ayo kita gambar histogram price dan log (price) untuk membandingkannya
4.1 Perbandingan menggunakan Log dan Tidak
Disini saya melakukan perbandingan dengan melihat elastisitasnya. Berikut koding nya:
#tanpa log
mod2 <- lm(price~sqft, data=br)
b1 <- coef(mod2)[[1]]
b2 <- coef(mod2)[[2]]
sqftx <- c(2000, 4000, 6000)
pricex <- b1+b2*sqftx^2
DpriceDsqft <- 2*b2*sqftx
short_e <- DpriceDsqft*sqftx/pricex
long_e <- short_e*(1/(1-b2))
short_e ; long_e## [1] 2.000328 2.000082 2.000036## [1] -0.02180257 -0.02179988 -0.02179939#dengan log
mod3 <- lm(log(price)~sqft, data=br)
b1 <- coef(mod3)[[1]]
b2 <- coef(mod3)[[2]]
sqftx <- c(2000, 4000, 6000)
pricex <- b1+b2*sqftx^2
DpriceDsqft <- 2*b2*sqftx
short_e <- DpriceDsqft*sqftx/pricex
long_e <- short_e*(1/(1-b2))
short_e ; long_e## [1] 1.986909 1.996711 1.998537## [1] 1.987727 1.997533 1.999359Jika diperhatikan, model regresi yang tanpa log menghasilkan elastisitas jangka pendek nya \(>1\) maka masih bisa dikatakan elastis atau nilai price lebih besar dibanding nilai sqft. Tetapi ketika elastisitasnya dalam jangka panjang menghasilkan nilainya \(<1\) yang mana merupakan inelastis atau nilai price lebih kecil dibanding nilai sqft. Dari hal ini, dapat disimpulkan bahwa ada kesalahan/ketidak akurat dalam perhitungannya. Sedangkan jika diperhatikan, model regresi yang menggunakan log menghasilkan elastisitas jangka pendek maupun jangka panjang nya \(>1\) maka masih masuk akal dan keduanya elasitis.
Seperti sebelumnya, kita tertarik dalam estimasi koefisien dan interpretasinya, pada nilai price yang sesuai, dan pada efek marginal dari kenaikan sqft pada price.
## (Intercept) sqft
## 1.083860e+01 4.112689e-04Koefisiennya adalah \(b_1 = 10,84\) dan \(b_2 = 0,00041\), menunjukkan bahwa peningkatan luas permukaan (sqft) sebuah apartemen sebesar satu unit (1 sqft) meningkatkan harga apartemen sebesar 0,041 persen. Jadi, untuk harga rumah sebesar $100.000, kenaikan 100 sqft akan menaikkan harga sekitar 100∗0,041 persen, yang sama dengan $4112,7. Secara umum, efek marjinal dari kenaikan \(X\) pada \(Y\) dalam Persamaan 5 adalah
\[ \begin{align} \tag{6} {dy\over dx}=\beta_2 * Y \end{align} \]
dan elastisitasnya adalah
\[ \begin{align} \tag{7} e= {dy\over dx} {X\over Y}=\beta_2 * X \end{align} \]
Kode berikutnya menunjukkan bagaimana menggambar kurva nilai pas dari model log-linear dan bagaimana menghitung efek marginal serta elastisitas harga median dalam dataset. Nilai yang pas dihitung menggunakan rumus
\[ \begin{align} \tag{8} \hat{Y}= e^{b_1+b_2*X} \end{align} \]
ordat <- br[order(br$sqft), ] #urutkan 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.9366934R memperbolehkan kita menghitung jumlah yang sama untuk beberapa pasang (sqft, price) sekaligus, seperti yang ditunjukkan dalam urutan berikut:
b1 <- coef(mod4)[[1]]
b2 <- coef(mod4)[[2]]
#ambil beberapa nilai untuk sqft:
sqftx <- c(2000, 3000, 4000)
#estimasi harga untuknya dan tambahkan satu lagi:
pricex <- c(100000, exp(b1+b2*sqftx))
#hitung ulang sqft untuk semua harga:
sqftx <- (log(pricex)-b1)/b2
#hitung dan cetak elastisitas:
(elasticities <- b2*sqftx) ## [1] 0.6743291 0.8225377 1.2338066 1.64507545 Variabel Indikator
Sebuah indikator, atau variabel binary menandai ada atau tidaknya beberapa atribut dari unit observasi, seperti jenis kelamin atau ras jika unit observasi adalah individu, atau lokasi jika unit observasi adalah sebuah rumah. Dalam dataset utown, variabelutown adalah 1 jika rumah dekat dengan universitas dan 0 jika sebaliknya. Berikut adalah model regresi linier sederhana yang melibatkan variabel utown:
\[ \begin{align} \tag{9} price_i= \beta_1+\beta_2*utwon_i \end{align} \]
Koefisien variabel semacam itu dalam model linier sederhana sama dengan selisih antara rata-rata harga dari dua kategori; koefisien intersep model pada Persamaan 9 sama dengan rata-rata harga rumah yang tidak dekat dengan universitas. Kita hitung dulu harga rata-rata untuk setiap kategori, yang dilambangkan dengan urutan kode price0bar dan price1bar:
data("utown", package="PoEdata")
price0bar <- mean(utown$price[which(utown$utown==0)])
price1bar <- mean(utown$price[which(utown$utown==1)])Hasilnya adalah: \(\bar{price}=\) 277.2416012 dekat ke universitas, dan \(\bar{price}=\) 215.7324948sebaliknya. Sekarang saya menunjukkan bahwa hasil yang sama menghasilkan koefisien model regresi dalam Persamaan 9:
## [1] 215.7325## [1] 61.50911Hasilnya adalah: \(\bar{price}=b_1=\) 215.7324948 untuk rumah non-universitas, dan \(\bar{price}=b_1+b_2=\) 277.2416012 untuk rumah universitas.
6 Simulasi Monte Carlo
Simulasi Monte Carlo menghasilkan nilai random untuk variabel dependen ketika koefisien regresi dan distribusi istilah random diberikan. Contoh berikut menentukan distribusi variabel independen dalam model biaya bahan makanan pada 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"))
Selanjutnya, kita hitung varians \(b_2\) dan plot fungsi kepadatannya.
\[ \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 yang sama \(b_1\), \(b_2\), dan error standard deviation, kita dapat menghasilkan satu set nilai untuk \(Y\), regresi \(Y\) pada \(X\), dan menghitung nilai estimasi untuk koefisien \(b_2\) dan standard error.
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)
seb2hat <- coef(mod6summary)[2,2]Hasilnya adalah \(b_2 = 11,64\) dan \(s_e (b_2) = 1,64\). Kekuatan dari simulasi Monte Carlo adalah kemungkinan pengulangan estimasi parameter regresi untuk sejumlah besar sampel yang dibuat secara otomatis. Jadi, kita dapat memperoleh sejumlah besar nilai untuk parameter, misalnya \(b_2\), dan kemudian menentukan karakteristik pengambilan sampelnya. Misalnya, jika mean dari nilai-nilai ini mendekati nilai yang diasumsikan semula \(b_2 = 10\), kita menyimpulkan bahwa estimator kita (metode estimasi parameter) tidak bias.
Kita akan menggunakan nilai \(X\) dalam dataset food, dan menghasilkan \(Y\) menggunakan model linier dengan \(b_1 = 100\) dan \(b_2 = 10\).
## Warning in data("food"): data set 'food' not foundN <- 40
sde <- 50
x <- food$income
nrsim <- 1000
b1 <- 100
b2 <- 10
vb2 <- numeric(nrsim)
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 \(b_2\), masing-masing adalah 9,974985 dan 1,152632. Gambar di bawah ini menunjukkan distribusi simulasi \(b_2\) dan teoritis.
plot(density(vb2))
curve(dnorm(x, mb2, seb2), col="red", add=TRUE)
legend("topright", legend=c("true", "simulated"),
lty=1, col=c("red", "black"))
