Email             :
Instagram     :
RPubs            :
Github           :
Telegram       :
Department  : [Business Statistics]
Address         :
                        


1 Introduction

Seperti yang telah saya sebutkan sebelumnya, Regression Models digunakan untuk menggambarkan hubungan antara variabel dengan memasang garis ke data yang diamati. Model regresi linear menggunakan garis lurus, sementara model regresi logistik dan nonlinear menggunakan garis melengkung. Regresi memungkinkan Anda untuk memperkirakan bagaimana variabel dependen berubah saat variabel independen berubah.

Di bagian ini kita akan memfokuskan cara menggunakan Simple Linear Regression (SLR) untuk memperkirakan hubungan antara dua variabel kuantitatif. Kita juga dapat menggunakan regresi linier sederhana untuk mengetahui:

  • Seberapa kuat hubungan antara dua variabel.
  • Nilai variabel dependen pada nilai tertentu dari variabel independen



Contoh:

Asumsikan Anda adalah peneliti sosial yang tertarik pada hubungan antara pendapatan dan kebahagiaan. Anda mensurvei 500 orang yang pendapatannya berkisar dari $15k sampai $75k dan minta mereka untuk mengurutkan kebahagiaan mereka berdasarkan skala 1 sampai 10. Variabel independen (pendapatan) dan variabel dependen (kebahagiaan) Anda berdua kuantitatif, sehingga Anda dapat melakukan analisis regresi untuk melihat apakah ada hubungan linier di antara mereka.

2 Model Umum

Model regresi linier sederhana mengasumsikan bahwa hubungan linier ada antara ekspektasi bersyarat dari dua variabel kuantitatif:

di mana,

  • \(Y\) dianggap sebagai variabel prediktor, penjelasan, atau dependen.
  • \(X\) dianggap sebagai respons, hasil, atau variabel independen.
  • \(\beta_1\) adalah parameter intercept atau koefisien
  • \(\beta_2\) adalah parameter kemiringan atau koefisien
  • \(\hat{Y}\) adalah estimasi nilai \(Y\) yang diberikan \(X\)
  • \(b_1\) dan \(b_2\) adalah nilai perkiraan koefisien

3 Regresi Linear Sederhana di 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 diinstall
library(devtools)                                    # muat librarynya
install_git("https://github.com/ccolonescu/PoEdata") # install `PoEdata`

Kemudian, pertimbangkan kumpulan data:

library(PoEdata)                                     # load the library PoEdata
library(DT)                                          # use for datatable fuction 
data("food", package="PoEdata")
datatable(food)

Selalu merupakan ide yang baik untuk secara visual memeriksa data dalam diagram sebar, yang dapat dibuat menggunakan fungsi plot(). Gambar di bawah adalah diagram sebar dari pengeluaran makanan terhadap pendapatan, menunjukkan 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")

About \(R^2\)

3.1 Memperkirakan Regresi Linear

Untuk data pengeluaran makanan, model regresi 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) yang mana digunakan dengan sendirinya tidak menunjukkan output apa pun; Berguna untuk memberi nama model, seperti mod1, then show the results using summary(mod1).

mod1 <- lm(food_exp ~ income, data = food)           # linear model
smod1<- summary(mod1)                                # get 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 hanya tertarik pada beberapa hasil regresi, seperti perkiraan koefisien, Anda dapat 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
b2
## [1] 10.20964

Dari hasilnya, kita dapat menulis regresi linier yang memperkirakan sebagai:

\[ \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 sedikit penting dalam model ekonometrik; kami sebagian besar tertarik pada parameter kemiringan, \(\beta_2\). Nilai perkiraan dari \(\beta_2\) menunjukkan bahwa pengeluaran makanan untuk keluarga rata-rata meningkat sebesar 10,209643 ketika pendapatan keluarga meningkat sebesar 1 unit, yang dalam hal ini $100. Fungsi R abline() menambahkan garis regfression line ke diagram sebar yang diplot sebelumnya:

plot(food$income, food$food_exp, 
     xlab="weekly income in $100", 
     ylab="weekly food expenditure in $", 
     type = "p")
abline(b1,b2)

Bagaimana seseorang dapat mengambil berbagai hasil regresi? Untuk mengambil hasil tertentu, Anda cukup merujuknya dengan nama objek, diikuti oleh tanda $ dan nama hasil yang ingin Anda ambil. Misalnya, jika kita ingin vektor koefisien dari mod1, kami menyebutnya sebagai mod1\$coefficients dan smod1\$coefficients:

mod1$coefficients
## (Intercept)      income 
##    83.41600    10.20964
smod1$coefficients
##             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
names(mod1)                                          # for more
names(smod1)                                         # for more

Seperti yang telah kita lihat sebelumnya, namun, beberapa hasil ini dapat diambil menggunakan fungsi tertentu, seperti coef(mod1), resid(mod1), fitted(mod1), and vcov(mod1).

3.2 Prediksi dengan Model SLR

Parameter regresi yang diperkirakan, \(b_1\) dan \(b_2\) memungkinkan kami untuk memprediksi pengeluaran makanan yang diharapkan untuk pendapatan tertentu. Yang perlu kita lakukan adalah menyambungkan 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 sebesar 100, karena data untuk pendapatan variabel adalah dalam ratusan dolar.)

\[ \begin{align} \tag{3} \hat{food_{exp}}&=83.416+10.20964*20\\ & = $ 287.608861 \end{align} \]

R, bagaimanapun, perhitungan ini bagi kita dengan fungsinya disebut predict(). Mari kita memperluas sedikit contoh untuk lebih dari satu pendapatan yang kita prediksi pengeluaran makanan, katakanlah pendapatan = $2000, $2500, dan $2700. Fungsi predict() dalam R mengharuskan nilai baru variabel independen diatur di bawah formulir tertentu, yang disebut bingkai data. Bahkan ketika kita hanya ingin memprediksi untuk satu pendapatan, kita membutuhkan struktur bingkai data yang sama. Dalam 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.0764

3.3 Menilai Koefisien Regresi

Koefisien regresi \(b_1\) dan \(b_2\) adalah variabel acak, karena tergantung pada sampel. Dalam hal ini, kami mengulangi sampel untuk menilai koefisien regresi. Mari kita membangun sejumlah subsample acak dari data makanan dan menghitung ulang \(b_1\) dan \(b_2\) . Subsample acak dapat dibangun menggunakan fungsi sample(), sebagai contoh berikut hanya mengilustrasikan untuk \(b_2\).

N <- nrow(food)   # returns the number of observations in the dataset
C <- 50           # desired number of subsamples
S <- 38           # desired sample size

sumb2 <- 0
for (i in 1:C){   # a loop over the number of subsamples
  set.seed(3*i)   # a different seed for each subsample  
  subsample <- food[sample(1:N, size=S, replace=TRUE), ]
  mod2 <- lm(food_exp~income, data=subsample)
  #sum b2 for all subsamples:
  sumb2 <- sumb2 + coef(mod2)[[2]]
}
print(sumb2/C, digits = 3)
## [1] 9.89

Hasilnya, \(b_2=9.88\), adalah rata-rata dari 50 estimasi \(b_2\).

3.4 Perkiraan Variansi dan Kovariansi

Banyak aplikasi memerlukan perkiraan varians dan kovarians koefisien regresi. R menyimpannya dalam matriks vcov():

(varb1 <- vcov(mod1)[1, 1])
## [1] 1884.442
(varb2 <- vcov(mod1)[2, 2])
## [1] 4.381752
(covb1b2 <- vcov(mod1)[1,2])
## [1] -85.90316

3.5 Tugas

  • Untuk hasil variansi dari b1, terlihat bahwa hasilnya sangat besar, yaitu 1884.442. Hal itu berarti terjadi penyimpangan yang sangat jauh dari rata-ratanya. Semakin besar nilai varians, semakin jauh data yang kita gunakan tersebar dari nilai rata-ratanya.

  • Untuk hasil variansi dari b2, terlihat bahwa hasilnya kecil, yaitu 4.381752. Hal itu berarti penyimpangan yang ternjadi hanya sedikit dari rata-ratanya. Semakin kecil nilai varians, semakin dekat data yang kita gunakan tersebar dari nilai rata-ratanya.

  • Untuk hasil kovariansi dari b1 dan b2, terlihat bahwa hasilnya adalah negatif, yaitu -85.90316. Hal itu berarti hubungan antarvariabelnya berlawanan, jika variabel yang satu naik, maka yang satu lagi turun, atau sebaliknya. ***

4 Hubungan Non-Linear

Kadang-kadang diagram plot sebar atau beberapa pertimbangan teoritis menunjukkan hubungan non-linear. Hubungan non-linear 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 dalam persamaan regresi dengan fungsi I(). Contoh berikut menggunakan dataset br dari paket PoEdata, yang mencakup harga jual dan atribut 1080 rumah di Baton Rouge, LA. harga adalah harga jual dalam dolar, dan sqft adalah area permukaan di 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)            # given values for sqft
pricex=b1+b2*sqftx^2                 # prices corresponding to given sqft 
DpriceDsqft <- 2*b2*sqftx            # marginal effect of sqft on price
elasticity=DpriceDsqft*sqftx/pricex 
b1; b2; DpriceDsqft; elasticity      # prints results
## [1] 55776.57
## [1] 0.0154213
## [1]  61.68521 123.37041 185.05562
## [1] 1.050303 1.631251 1.817408

Kami ingin sekarang menggambar diagram sebar dan melihat bagaimana fungsi kuadrat sesuai dengan data. Potongan kode berikutnya menyediakan dua alternatif untuk membangun grafik seperti itu. Yang pertama hanya menggambar fungsi kuadrat pada diagram sebar, menggunakan fungsi R curve(); kedua menggunakan garis fungsi, yang mengharuskan urutan kumpulan data dalam meningkatkan nilai sqft sebelum model regresi dievaluasi, sedemikian rupa sehingga nilai yang sesuai 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 dipasang:

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")

Model log-linear meregresi log variabel dependen pada ekspresi linear dari variabel independen (kecuali ditentukan lain, notasi log adalah singkatan dari logaritma alami, mengikuti konvensi biasa di bidang ekonomi):

\[ \begin{align} \tag{5} Log(Y)=\beta_1+\beta_2*X+e\\ \end{align} \]

Salah satu alasan untuk menggunakan log variabel independen adalah untuk membuat distribusinya lebih dekat ke distribusi normal. Mari kita menggambar histogram harga dan log (harga) untuk membandingkannya

hist(br$price, col='grey')

hist(log(br$price), col='grey')

Kami tertarik, seperti sebelumnya, dalam perkiraan koefisien dan interpretasinya, dalam nilai harga yang sesuai, dan dalam efek marginal dari peningkatan sqft pada harga.

mod3 <- lm(log(price)~sqft, data=br)
mod3$coefficients
##  (Intercept)         sqft 
## 1.083860e+01 4.112689e-04

Koefisien \(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. Dengan demikian, untuk harga $100,000, meningkat 100 sqft akan meningkatkan harga sekitar 100∗0,041 persen, yang sama dengan $4112.7. Secara umum, efek marjinal dari peningkatan \(X\) on \(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 untuk harga rata-rata dalam set data. Nilai yang dipasang 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")

pricex<- median(br$price)
sqftx <- (log(pricex)-coef(mod4)[[1]])/coef(mod4)[[2]]
(DyDx <- pricex*coef(mod4)[[2]])
## [1] 53.46495
(elasticity <- sqftx*coef(mod4)[[2]])
## [1] 0.9366934

R memungkinkan kita untuk menghitung jumlah yang sama untuk beberapa pasangan (sqft, harga) pada satu waktu, 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.6450754

5 Indicator Variables

Indikator, atau variabel biner menandai keberadaan atau tidak adanya beberapa atribut unit observasional, seperti jenis kelamin atau ras jika unit observasional adalah individu, atau lokasi jika unit observasional adalah rumah. Pada dataset utown, variabel utown adalah 1 jika rumah dekat dengan universitas dan 0 sebaliknya. Berikut adalah model regresi linear sederhana yang melibatkan variabel u town:

\[ \begin{align} \tag{9} price_i= \beta_1+\beta_2*utwon_i \end{align} \]

The coefficient of such a variable in a simple linear model is equal to the difference between the average prices of the two categories; the intercept coefficient of the model in Equation 9 is equal to the average price of the houses that are not close to university. Let us first calculate the average prices for each category, wich are denoted in the following sequence of code price0bar and price1bar:

data(utown)
price0bar <- mean(utown$price[which(utown$utown==0)])
price1bar <- mean(utown$price[which(utown$utown==1)])

Hasilnya: \(\bar{price}=\) 277.2416012 dekat dengan universitas, dan \(\bar{price}=\) 215.7324948bagi mereka yang tidak dekat. Saya sekarang menunjukkan bahwa hasil yang sama menghasilkan koefisien model regresi dalam Persamaan 9 :

mod4 <- lm(price~utown, data=utown)
b1 <- coef(mod4)[[1]] 
b2 <- coef(mod4)[[2]]
b1
## [1] 215.7325
b2
## [1] 61.50911

Hasilnya: \(\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 acak untuk variabel dependen ketika koefisien regresi dan distribusi istilah acak diberikan. Contoh berikut berusaha untuk menentukan distribusi variabel independen dalam 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"))

Selanjutnya, kita menghitung varian \(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)

Now, with the same values of \(b_1\) , \(b_2\) , and error standard deviation, we can generate a set of values for \(Y\) , regress \(Y\) on \(X\) , and calculate an estimated values for the coefficient \(b_2\) and its 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) #the summary contains the standard errors
seb2hat <- coef(mod6summary)[2,2]

Hasilnya \(b_2=11.64\) and \(s_e(b_2)=1.64\). Kekuatan simulasi Monte Carlo adalah, bagaimanapun, kemungkinan mengulangi estimasi parameter regresi untuk sejumlah besar sampel yang dihasilkan secara otomatis. Dengan demikian, kita dapat memperoleh sejumlah besar nilai untuk parameter, katakanlah \(b_2\) , dan kemudian menentukan karakteristik pengambilan sampelnya. Misalnya, jika nilai rata-rata ini mendekati nilai yang awalnya tersederakan \(b_2= 10\), kami menyimpulkan bahwa estimator kami (metode memperkirakan parameter) tidak bias.

Kita akan menggunakan kali ini nilai-nilai \(X\) dalam kumpulan data makanan, 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) #stores the estimates of 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)

Simpangan rata-rata dan baku dari perkiraan 40 nilai \(b_2\) masing-masing adalah 9,974985 dan 1,152632. Gambar bellow menunjukkan simulasi distribusi \(b_2\) dan yang 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"))

hist(vb2, prob=TRUE, ylim=c(0,.4))
curve(dnorm(x, mean=mb2, sd=seb2), col="red", add=TRUE)

rm(list=ls()) # Caution: this clears the Environment