Econometrics
Simple Linear Regression
Vanessa Supit
February 09, 2021
Email : vanessasupit0910@gmail.com
RPubs : https://rpubs.com/vanessasupit/
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 disebutkan sebelumnya, Model Regresi digunakan untuk memperjelas hubungan antara variabel dengan menyesuaikan garis ke data yang diamati. Model regresi linear menggunakan garis lurus, sedangkan Model regresi logistik dan nonlinear menggunakan garis kurva. Regresi membuat anda dapat mengestimasi bagaimana variabel dependen berubah saat variabel independentnya berubah.
Pada bagian ini kita akn fokus cara menggunakan Simple Linear Regresi (SLR) atau Regresi Linear Sederhana untuk mengestimasi hubungan antara dua kuantiatif variabel. kita dapat juga mengunakan SLR untuk mengetahui :
- Kuat pengaruh hubungan antar variabel
- Nilai variabel dependen saat nilai variabel independent tertentu
Contoh:
Asumsikan anda adalah peneliti sosial yang tertarik dalam hubungan antara pendapatan dan kebahagiaan. Ansa mensurvei 500 orang yang rentang gajinya dari $15k sampai $75k dan meminta mereka untuk memberi nilai kebahagiaan mereka dari 1 sampai 10. Variabel bebas(Gaji) dan variabel terikat (bahagia) mu adalah kuantitatif, sehingga anda dapat melakukan analisa regresi untuk melihat apakah ada hubungan linier antara mereka.
2 Model Umum
Model Regresi Linear Sederhana mengasumsikan bahwa terdapat hubungan antara kondisi harapatan dari dua vairabel kuantitatif :
dimana,
- \(Y\) melambangkan respon, hasil, atau variabel terikat.
- \(X\) melambangkan prediktor, penjelas, atau variabel bebas.
- \(\beta_1\) adalah parameter atau koefisien perpotongan
- \(\beta_2\) adalah parameter atau koefisien gradien
- \(\hat{Y}\) adalah nilai estimasi dari \(Y\) untuk \(X\)
- \(b_1\) dan \(b_2\) adalah nilai koefisien estimasi
3 Regresi Linear Sederhana dalam R
Data untuk contoh ini disimpan dalam package R PoEdata. Pertama, anda harus menginstall package “PoEdata”, ketik script berikut di konsol R studio :
install.packages("devtools") # jika belum diinstal
library(devtools) # muat library
install_git("https://github.com/ccolonescu/PoEdata") # install `PoEdata`Lalu, pertimbangkan kumpulan data:
library(PoEdata) # muat library PoEdata
library(DT) # gunakan untuk fungsi datatable
data("food", package="PoEdata")
datatable(food)Merupakan ide yang bagus untuk memvisualisasikan data dalam diagram scatter, dimana dapat dibuat menggunaan fungsi plot(). gambar berikutan merupakan diagram scatter dari pengeluaran makanan terhadap pendapatan, anggaplah 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 Mengestimasikan Regresi Linear
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 mengestimasi model Regresi Linear lm(y~x, data) dimana, digunakan dirinya sendiri tidak menunjukkan output apapun; Sangat membantu jikak memberikan namauntuk model, seperti mod1, lalu menunjukkan hasilnya menggunakan summary(mod1).
mod1 <- lm(food_exp ~ income, data = food) # Model Linear
smod1<- summary(mod1) # mendapat kesimpulan
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-05Kalau kamu hanya tertarik pada beberapa hasil regresi, seperti koefisien estimasi, anda bisa mendapatkannya menggunakan fungsi tertentu, seperti fungsi coef().
## [1] 83.416## [1] 10.20964Dari hasil ini, kita dapat menulis estimasi linier regresi sebagai berikut:
\[ \begin{align} \tag{2} \hat{Y}&=b_1+b_2*X \\ \hat{food_{exp}}&=83.416+10.20964*income \end{align} \]
Parameter perpotongan, \(\beta_1\) , biasanya tidak terlalu penting dalam model ekometrika; kita lebih tertarik dalam parameter gradien, \(\beta_2\). Nilai estimasi dari \(\beta_2\) menunjukkan bahwa pengeluaran makanan untuk rata-rata keluarga meningkat 10.209643 disaat pendapatan keluarga meningkat 1 unit, yang mana dalam kasus ini $100. Fungsi R abline() menambahkan garis regresi ke dalam diagram scatter plot sebelumnya di bawah 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 hasil regresi? Untuk mengambil hasil tertentu Anda cukup merujuknya dengan nama objek, diikuti dengan simbol $ dan nama hasil yang ingin Anda ambil. Misalnya, jika kita ingin vektor koefisien dari mod1, kami menyebutnya sebagai mod1\$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-05Namun, seperti yang telah kita lihat sebelumnya, beberapa dari hasil ini dapat diambil menggunakan fungsi tertentu, seperti coef(mod1), resid(mod1), fitted(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 penghasilan $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, Namun, apakah perhitungan ini untuk kita dengan fungsinya dipanggil predict(). Mari kita sedikit memperluas contoh ke lebih dari satu pendapatan yang kita perkirakan pengeluaran makanannya, katakanlah pendapatan = $2000, $2500, dan $2700. Fungsi predict() di R mensyaratkan bahwa 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, sekumpulan 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 # menampilkan hasil## 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 yang diilustrasikan contoh berikut hanya untuk \(b_2\).
N <- nrow (food) # mengembalikan jumlah observasi dalam dataset
C <- 50 # jumlah subsampel yang diinginkan
S <- 38 # ukuran sampel yang diinginkan
sumb2 <- 0
for (i in 1:C){ # loop di atas jumlah subsampel
set.seed (3 * i) # seed berbeda untuk setiap subsampel
subsample <- food[sample(1:N, size=S, replace=TRUE), ]
mod2 <- lm(food_exp~income, data=subsample)
#menjumlahkan b2 untuk semua subsampel:
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 [Tugas 1] Perkiraan Variansi dan Kovariansi
Banyak aplikasian membutuhkan estimasi dari variansi dan kovariansi koefisien regresi. R menempatkan mereka 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}
\] Dalam R, variabel independen yang melibatkan operator matematika dapat dimasukkan ke dalam persamaan regresi dengan fungsi I(). Contoh berikut menggunakan dataset br dari package PoEdata, yang mencakup harga jual dan atribut 1.080 rumah di Baton Rouge, LA. harga adalah harga jual dalam dolar, dan sqft adalah luas permukaan dalam sqft.
data(br)
mod2 <- lm(price~I(sqft^2), data=br)
b1 <- coef(mod2)[[1]]
b2 <- coef(mod2)[[2]]
sqftx=c(2000, 4000, 6000) # nilai yg 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 # menampilkan hasil## [1] 55776.57## [1] 0.0154213## [1] 61.68521 123.37041 185.05562## [1] 1.050303 1.631251 1.817408Kami ingin menggambar diagram sebar 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 membutuhkan pengurutan kumpulan data dalam meningkatkan 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")
#add the quadratic curve to the scatter plot:
curve(b1+b2*x^2, col="red", add=TRUE)
Cara alternatif untuk menggambar kurva yang pas:
ordat <- br[order(br$sqft), ] #sorts the dataset after `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(mod2)~ordat$sqft, col="red")4.1 [Tugas 2] Model log-linear
Model log-linear meregresi log dari variabel dependen pada ekpresi linier dari variabel bebas (kecuali ditentukan lain, 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 variabel bebas adalah agar distribusinya mendekati distribusi normal. Mari kita menggambar histogram harga dan log (harga) untuk membandingkannya
seperti sebelumnya kita tertarik 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-04Koefesiennya adalah \(b_1= 10.84\) dan \(b_2= 0.00041\), menunjukkan bahwa peningkatan luas permukaan (sqft) apartemen sebesar satu unit (1 sqft) meningkatkan harga apartemen sebesar 0,041 persen. Jadi, untuk harga rumah sebesar $100,000, peningkatan 100 sqft akan meningkatkan harga sekitar 100 ∗ 0,041 persen, yang sama dengan $4112.7. Secara umum, efek marginal dari peningkatan \(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} \]
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), ] #order the 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")
abline(lm(ordat$price~ordat$sqft))
text("hahaha")## Warning in xy.coords(x, y, recycle = TRUE, setLab = FALSE): NAs introduced by
## coercion
Bisa kita lihat dari perbandingan di atas bahwa grafik logaritma (warna biru) cocok dengan data-data yang tersedia, sedangkan grafik regresi sederhana(warna hitam) hanya cocok dengan data-data yang bernilai kecil namun ketika datanya mulai menyebar atau meluas, pencilannya semakin banyak, prediksi yang kita akan menyimpang. pada grafik logaritma lebih memungkinkan karena pencilannya kecil karena bentuk dari grafik logaritma lebih bersesuaian 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]]
#pick a few values for sqft:
sqftx <- c(2000, 3000, 4000)
#estimate prices for those and add one more:
pricex <- c(100000, exp(b1+b2*sqftx))
#re-calculate sqft for all prices:
sqftx <- (log(pricex)-b1)/b2
#calculate and print elasticities:
(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, variabel utown bernilai 1 jika sebuah rumah dekat dengan universitas dan 0 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 dari dua 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 dalam urutan kode berikut price0bar dan price1bar:
data(utown)
price0bar <- mean(utown$price[which(utown$utown==0)])
price1bar <- mean(utown$price[which(utown$utown==1)])Hasilnya adalah: \(\bar{price}=\) 277.2416012 mendekati university, dan \(\bar{price}=\) 215.7324948bagi mereka yang tidak dekat. 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-university, dan \(\bar{price}=b_1+b_2=\) 277.2416012 untuk rumah university.
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 kita menghitung variansi dari \(b_2\) dan plot 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 dari \(b_1\) , \(b_2\) , dan standar deviasi kesalahan, kita dapat menghasilkan satu set nilai untuk \(Y\) , regresi \(Y\) pada \(X\) , dan menghitung nilai perkiraan untuk koefisien \(b_2\) dan standar kesalahannya
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) #the summary contains the standard errors
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 bisa mendapatkan sejumlah besar nilai untuk suatu parameter, misalnya \(b_2\) , dan kemudian menentukan karakteristik pengambilan sampelnya. Misalnya, jika mean dari nilai-nilai ini dekat dengan nilai yang diasumsikan semula \(b_2= 10\), kami menyimpulkan bahwa estimator kami (metode estimasi parameter) tidak bias.
Kami akan menggunakan kali ini nilai dari \(X\) dalam kumpulan data 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 estimasi 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 dari \(b_2\) masing-masing adalah 9.974985 dan 1.152632. Gambar di bawah ini menunjukkan distribusi simulasi \(b_2\) dan 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"))
