Econometrics
Simple Linear Regression
Jocelyn Irene Gani
March 02, 2021
Email : jirene113@gmail.com
Instagram : https://www.instagram.com/irenegani
RPubs : https://rpubs.com/irenegani/
Department : Business Statistics
Address : ARA Center, Matana University Tower
Jl. CBD Barat Kav, RT.1, Curug Sangereng, Kelapa Dua, Tangerang, Banten 15810.
1 Perkenalan
Seperti yang telah saya sebutkan sebelumnya, * Model Regresi * digunakan untuk menggambarkan hubungan antara variabel dengan menyesuaikan garis ke data yang diamati. Model regresi linier menggunakan garis lurus, sedangkan model regresi logistik dan nonlinier menggunakan garis lengkung. Regresi memungkinkan Anda memperkirakan bagaimana variabel dependen berubah saat variabel independen berubah.
Pada bagian ini kita akan fokus pada bagaimana menggunakan ** Regresi Linier Sederhana (SLR) ** untuk memperkirakan hubungan antara dua variabel kuantitatif. Kita juga bisa menggunakan regresi linier sederhana untuk mengetahui:
- Seberapa kuat hubungan antara dua variabel.
- Nilai variabel terikat pada nilai tertentu dari variabel bebas
Contoh:
Asumsikan Anda adalah peneliti sosial yang tertarik pada hubungan antara pendapatan dan kebahagiaan. Anda mensurvei 500 orang yang pendapatannya berkisar dari $ 15.000 hingga $ 75.000 dan meminta mereka untuk mengurutkan kebahagiaan mereka dalam skala dari 1 hingga 10. Variabel bebas (pendapatan) dan variabel terikat (kebahagiaan) Anda sama-sama bersifat kuantitatif, sehingga Anda dapat melakukan analisis regresi untuk melihat apakah terdapat hubungan linier di antara keduanya.
2 Model Umum
Model regresi linier sederhana mengasumsikan bahwa ada hubungan linier antara ekspektasi bersyarat dari dua variabel kuantitatif:
dimana,
- \(Y\) dianggap sebagai prediktor, penjelas, atau variabel dependen.
- \(X\) dianggap sebagai respon, hasil, atau variabel independen.
- \(\beta_1\) adalah parameter atau koefisien intersep
- \(\beta_2\) adalah parameter atau koefisien kemiringan
- \(\hat{Y}\) adalah perkiraan nilai \(Y\) diberikan \(X\)
- \(b_1\) dan \(b_2\) adalah nilai estimasi koefisien
3 Regresi Linier Sederhana dalam R
Data untuk contoh ini disimpan dalam paket R PoEdata. Pertama, Anda harus menginstal paket “PoEdata”, ketik baris skrip berikut di jendela Konsol RStudio:
# install.packages("devtools") # jika belum terpasang
library(devtools) # memuat library
install_git("https://github.com/ccolonescu/PoEdata") # install `PoEdata`Kemudian, pertimbangkan kumpulan data:
library(PoEdata) # memuat library PoEdata
library(DT) # digunakan untuk fungsi dataable
data("food", package="PoEdata")
datatable(food)Itu selalu merupakan ide yang baik untuk memeriksa data secara visual dalam diagram pencar, yang dapat dibuat menggunakan fungsi plot (). Gambar di bawah ini adalah diagram sebar pengeluaran makanan terhadap pendapatan, yang menunjukkan bahwa terdapat 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")
3.1 Memperkirakan Regresi Linear
Untuk data pengeluaran pangan, akan dibuat model regresi
\[ \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 linier adalah lm (y ~ x, data) yang, digunakan dengan sendirinya tidak menunjukkan keluaran apa pun; Berguna untuk memberi nama pada model, seperti mod1, lalu tampilkan hasil menggunakanringkasan (mod1).
mod1 <- lm(food_exp ~ income, data = food) # model linier
smod1<- summary(mod1) # dapatkan ringkasan
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 Anda tertarik hanya pada beberapa hasil regresi, seperti koefisien yang diperkirakan, Anda bisa mengambilnya menggunakan fungsi tertentu, seperti fungsi coef ().
b1 <- coef(mod1)[[1]] # the estimated value of β1
b2 <- coef(mod1)[[2]] # the estimated value of β2
b1## [1] 83.416## [1] 10.20964Dari hasil tersebut, kita dapat menuliskan estimasi regresi linier kita sebagai:
\[ \begin{align} \tag{2} \hat{Y}&=b_1+b_2*X \\ \hat{food_{exp}}&=83.416+10.20964*income \end{align} \]
Parameter intersep, $ beta_1 $, biasanya tidak terlalu penting dalam model ekonometrik; kami sangat tertarik dengan parameter kemiringan, $ beta_2 $. Perkiraan nilai $ beta_2 $ menunjukkan bahwa pengeluaran makanan untuk rata-rata keluarga meningkat sebesar 10.209643 ketika pendapatan keluarga meningkat sebesar 1 unit, yang dalam hal ini adalah $ 100. Fungsi R abline () menambahkan garis regfression ke diagram sebar yang sebelumnya diplot di bawah:
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 Anda cukup merujuknya dengan nama objek, diikuti dengan tanda $ dan nama hasil yang ingin Anda 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 telah kita lihat sebelumnya, bagaimanapun, beberapa dari hasil ini dapat diambil menggunakan fungsi tertentu, seperti coef (mod1), resid (mod1), pas (mod1), dan vcov (mod1).
3.2 Prediksi dengan Model SLR
Parameter regresi yang diperkirakan, $ b_1 $ dan $ b_2 $ memungkinkan kita untuk memprediksi pengeluaran makanan yang diharapkan untuk pendapatan tertentu. Yang perlu kita lakukan adalah memasukkan nilai parameter yang diperkirakan 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} \]
R, bagaimanapun, melakukan perhitungan ini untuk kita dengan fungsinya yang disebut predict (). Mari kita tambahkan sedikit contoh ke lebih dari satu pendapatan yang kita prediksi pengeluaran makanannya, katakanlah pendapatan = $ 2000, $ 2500, dan $ 2700. Fungsi predict () di R mengharuskan nilai baru dari variabel independen diatur dalam bentuk tertentu, yang disebut kerangka data. Bahkan ketika kita hanya ingin memprediksi untuk satu pendapatan, kita membutuhkan struktur kerangka data yang sama. Di R, satu set angka disatukan menggunakan struktur c (). Urutan berikut menunjukkan contoh ini.
mod1 <- lm(food_exp~income, data=food)
newx <- data.frame(income = c(20, 25, 27))
yhat <- predict(mod1, newx)
yhat # prints the result## 1 2 3
## 287.6089 338.6571 359.07643.3 Nilai Koefisien Regresi
Koefisien regresi $ b_1 $ dan $ b_2 $ adalah variabel acak, karena bergantung pada sampel. Dalam kasus ini, kami mengulang sampel untuk menilai koefisien regresi. Mari kita buat sejumlah subsampel acak dari data makanan dan hitung ulang $ b_1 $ dan $ b_2 $. Sebuah subsampel acak dapat dibangun menggunakan fungsi sample (), seperti contoh berikut hanya menggambarkan untuk $ b_2 $.
N <- nrow(food) # mengembalikan jumlah observasi dalam kumpulan data
C <- 50 # jumlah subsampel yang diinginkan
S <- 38 # ukuran sampel yang diinginkan
sumb2 <- 0
for (i in 1:C){ # satu putaran di atas jumlah subsampel
set.seed(3*i) # benih yang berbeda untuk setiap subsampel
subsample <- food[sample(1:N, size=S, replace=TRUE), ]
mod2 <- lm(food_exp~income, data=subsample)
#jumlah b2 untuk semua sub-sampel:
sumb2 <- sumb2 + coef(mod2)[[2]]
}
print(sumb2/C, digits = 3)## [1] 9.89Hasilnya, \(b_2=9.88\), adalah rata-rata dari 50 perkiraan \(b_2\).
3.4 Perkiraan Var dan Cov
Banyak aplikasi memerlukan perkiraan varians dan kovarians dari koefisien regresi. R menyimpannya dalam matriks vcov ():
## [1] 1884.442paste("Variansi dari b1 tinggi yang mengindikasikan bahwa titik data tersebar disekitar rerata (nilai ekspetasi) dan dari satu sama lainnya")## [1] "Variansi dari b1 tinggi yang mengindikasikan bahwa titik data tersebar disekitar rerata (nilai ekspetasi) dan dari satu sama lainnya"## [1] 4.381752paste("Variansi dari b2 kecil yang mengindikasikan bahwa titik data condong sangat dekat dengan nilai rerata (nilai ekspetasi) dan dari satu sama lainnya")## [1] "Variansi dari b2 kecil yang mengindikasikan bahwa titik data condong sangat dekat dengan nilai rerata (nilai ekspetasi) dan dari satu sama lainnya"## [1] -85.90316paste("Kovariansi memiliki nilai negatif yang berarti terdapat hubungan namun terbalik dimana ketika X membesar maka Y mengecil ")## [1] "Kovariansi memiliki nilai negatif yang berarti terdapat hubungan namun terbalik dimana ketika X membesar maka Y mengecil "4 Hubungan Non-Linear
Terkadang diagram diagram sebar atau beberapa pertimbangan teoretis 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}
\] Pada R, variabel bebas yang melibatkan operator matematika dapat dimasukkan ke dalam persamaan regresi dengan fungsi I (). Contoh berikut menggunakan kumpulan data br dari paket PoEdata, yang mencakup harga jual dan atribut 1080 rumah di Baton Rouge, LA. harga adalah harga jual dalam dolar, dan kaki persegi adalah luas permukaan dalam kaki persegi.
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 sqft
pricex=b1+b2*sqftx^2 # harga yang sesuai dengan sqft yang diberikan
DpriceDsqft <- 2*b2*sqftx # efek marjinal sqft terhadap harga
elasticity=DpriceDsqft*sqftx/pricex
b1; b2; DpriceDsqft; elasticity # hasil cetakan## [1] 55776.57## [1] 0.0154213## [1] 61.68521 123.37041 185.05562## [1] 1.050303 1.631251 1.817408Kami ingin menggambar diagram pencar dan melihat bagaimana fungsi kuadrat cocok dengan data. Potongan kode berikutnya memberikan dua alternatif untuk membuat grafik seperti itu. Yang pertama hanya menggambar fungsi kuadrat pada diagram pencar, menggunakan fungsi R 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")
#tambahkan kurva kuadrat ke plot pencar:
curve(b1+b2*x^2, col="red", add=TRUE) 
Cara alternatif untuk menggambar kurva yang pas:
ordat <- br[order(br$sqft), ] #mengurutkan set data 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")4.1 Model log-linear meregresi log dari variabel dependen pada ekspresi linier dari variabel independen (kecuali ditentukan lain, notasi log adalah singkatan dari natural logarithm, 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 variabel bebas adalah agar distribusinya mendekati distribusi normal. Mari kita menggambar histogram harga dan log (harga) untuk membandingkannya
Kami tertarik, seperti sebelumnya, pada estimasi koefisien dan interpretasinya, pada nilai harga yang sesuai, dan pada efek marginal dari peningkatan sqft pada harga.
## (Intercept) sqft
## 1.083860e+01 4.112689e-04Koefisiennya adalah $ b_1 = 10,84 $ dan $ b_2 = 0,00041 \(, menunjukkan bahwa peningkatan luas permukaan (kaki persegi) sebuah apartemen sebesar satu unit (1 kaki persegi) meningkatkan harga apartemen sebesar 0,041 persen. Jadi, untuk harga rumah sebesar \ \) 100.000, kenaikan 100 kaki persegi akan menaikkan harga sekitar 100 ∗ 0,041 persen, yang sama dengan $ 4112,7. Secara umum, efek marginal dari kenaikan $ X $ pada $ Y $ dalam Persamaan 5 adalah
\[ \begin{align} \tag{6} {dy\over dx}=\beta_2 * Y \end{align} \]
dan elastisitasnya
\[ \begin{align} \tag{7} e= {dy\over dx} {X\over Y}=\beta_2 * X \end{align} \]
Baris kode berikutnya menunjukkan cara menggambar kurva nilai yang dipasang dari model loglinear dan cara menghitung efek marginal dan elastisitas harga median dalam kumpulan data. Nilai yang sesuai di sini dihitung menggunakan rumus
\[ \begin{align} \tag{8} \hat{Y}= e^{b_1+b_2*X} \end{align} \]
ordat <- br[order(br$sqft), ] #memesan set data
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")
Bisa dilihat dari perbandingan di atas bahwa grafik logaritma (warna biru) cocok dengan data yang tersedia, sedangkan grafik regresi sederhana (warna hitam) hanya cocok dengan data-data yang bernilai kecil namun ketika datanya mulai menyebar, maka pencilannya semakin banyak, prediksi yang kita punya akan menyimpang. Pada grafik logaritma lebih memungkinkan karena pencilannya kecil karena bentuk dari grafik logaritma lebih sesuai dengan penyebaran data.
pricex<- median(br$price)
sqftx <- (log(pricex)-coef(mod4)[[1]])/coef(mod4)[[2]]
(DyDx <- pricex*coef(mod4)[[2]])## [1] 53.46495## [1] 0.9366934R memungkinkan kita menghitung jumlah yang sama untuk beberapa pasang (sqft, harga) sekaligus, seperti yang ditunjukkan dalam urutan berikut:
b1 <- coef(mod4)[[1]]
b2 <- coef(mod4)[[2]]
#pilih beberapa nilai untuk sqft:
sqftx <- c(2000, 3000, 4000)
#perkirakan harga untuk itu 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 biner 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 harga rata-rata kedua kategori; koefisien intersep model pada Persamaan 9 sama dengan harga rata-rata rumah yang tidak dekat dengan universitas. Mari kita hitung dulu harga rata-rata untuk setiap kategori, yang dilambangkan dengan urutan kode price0bar danprice1bar berikut:
data(utown)
price0bar <- mean(utown$price[which(utown$utown==0)])
price1bar <- mean(utown$price[which(utown$utown==1)])Hasilnya adalah: $ bar {price} = $ 277.2416012 dekat dengan universitas, dan $ bar {price} = $r price0baruntuk yang tidak tutup. 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 = $r b1 + b2 untuk rumah universitas.
6 Simulasi Monte Carlo
Simulasi Monte Carlo menghasilkan nilai acak untuk variabel dependen ketika koefisien regresi dan distribusi istilah acak diberikan. Contoh berikut berupaya menentukan distribusi variabel independen dalam model pengeluaran 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, kami menghitung varian $ b_2 $ dan memplot fungsi kepadatan yang sesuai.
\[ \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 deviasi standar kesalahan, kita dapat menghasilkan satu set nilai untuk $ Y $, regresi $ Y $ pada $ X $, dan menghitung nilai perkiraan untuk koefisien $ b_2 $ dan kesalahan standarnya.
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) #ringkasan berisi kesalahan standar
seb2hat <- coef(mod6summary)[2,2]Hasilnya adalah $ b_2 = 11,64 $ dan $ s_e (b_2) = 1,64 $. Kekuatan dari simulasi Monte Carlo adalah, bagaimanapun, 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 $, kami menyimpulkan bahwa estimator kami (metode estimasi parameter) tidak bias.
Kita akan menggunakan kali ini nilai $ X $ dalam dataset makanan, dan menghasilkan $ Y $ menggunakan model linier 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 perkiraan 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 deviasi standar dari perkiraan 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"))
