Email         : juliansalomo2@gmail.com
RPubs       : https://rpubs.com/juliansalomo/
Department : Business Statistics
Address : ARA Center, Matana University Tower
   Jl. CBD Barat Kav, RT.1, Curug Sangereng, Kelapa Dua, Tangerang, Banten 15810.

pacman::p_load(tidyr,
               dplyr,
               tidyverse,
               stats,
               car,
               DT,
               broom,
               ggplot2,
               PoEdata,
               effects,
               knitr
               )

1 Regresi Linier Sederhana VS Regresi Linier Berganda

Apa perbedaan regressi Linear Sederhana dan Berganda, jelaskan dengan contoh!

Regresi Linier Sederhana hanya memiliki satu variabel terikat dan bebas sedangkan Regresi Linier berganda memiliki 1 variabel terikat dan lebih dari 1 variabel bebas.

Contoh Regresi Linier Sederhana adalah sebagai berikut:

data("food")
datatable(food,
            extensions = 'FixedColumns',
            option = list(scrollX = TRUE, fixedColumns = TRUE)
          )

Data ini yang akan kita regresikan dengan fungsi lm(y~x).

linmod <- lm(food_exp ~ income, data = food)
summary(linmod)

## 
## 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

Dari bagian koefisien bisa kita lihat bahwa hanya ada satu variabel bebas yaitu variabel income

Selanjutnya contoh untuk Regresi Linier Berganda adalah sebagai berikut:

data("andy")
datatable(andy,
            extensions = 'FixedColumns',
            option = list(scrollX = TRUE, fixedColumns = TRUE)
          )

Data ini yang akan kita regresikan dengan fungsi lm(y~x1+x2).

linmod <- lm(sales ~ price + advert, data = andy)
summary(linmod)

## 
## Call:
## lm(formula = sales ~ price + advert, data = andy)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -13.4825  -3.1434  -0.3456   2.8754  11.3049 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept) 118.9136     6.3516  18.722  < 2e-16 ***
## price        -7.9079     1.0960  -7.215 4.42e-10 ***
## advert        1.8626     0.6832   2.726  0.00804 ** 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 4.886 on 72 degrees of freedom
## Multiple R-squared:  0.4483, Adjusted R-squared:  0.4329 
## F-statistic: 29.25 on 2 and 72 DF,  p-value: 5.041e-10

Dari bagian koefisien bisa kita lihat bahwa ada lebih dari satu variabel bebas, disini saya menggunakan dua variabel, yaitu variabel price dan advert

2 Analisis Regresi Linier Sederhana

Lakukan analisis regresi linear sederhana dalam ilmu ekonometrik!

Dalam analisis regresi linier sederhana, saya akan menggunakan data insur dari package PoEdata. Berikut data yang saya gunakan

data(insur)
datatable(insur,
          caption = htmltools::tags$caption(
            style = 'caption-side: bottom; text-align: center;',
            htmltools::em('Table Households Insurance and Income.')),
            extensions = 'FixedColumns',
            option = list(scrollX = TRUE, fixedColumns = TRUE)
          )

Data ini merupakan data dari hasil pengamatan untuk 20 rumah tangga. Yang diamati adalah pendapatan dari tiap rumah tangga dan biaya yang digunakan untuk insuransi. Untuk tahap awal kita dapat memulai dengan mengecek plot dari data kita. Kita bisa menggunakan fungsi plot()

plot(insur$income, insur$insurance,
     ylim = c(0, max(insur$insurance)),
     xlim = c(0, max(insur$income)),
     ylab = "Insurance in $1000",
     xlab = "Income in $1000"
     )

2.1 Mengestimasi Regresi Linier

Pertama-tama kita mulai dari mengestimasi model linier kita. Dalam data insur, yang menjadi variabel terikatnya adalah insurance sehingga model variabel kita adalah sebagai berikut

\[ \begin{align} Y &= \beta_1 + \beta_2 * X + e \\ Insurance &= \beta_1 + \beta_2 * income + e \end{align} \]

Selanjutnya kita dapat mencari koefisiennya dengan memasukkan model linier kita lalu mencari summary().

lm <- lm(insurance ~ income, data = insur)
slm <- summary(lm)
slm

## 
## Call:
## lm(formula = insurance ~ income, data = insur)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -24.228 -10.766   2.456  11.295  21.739 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)   6.8550     7.3835   0.928    0.365    
## income        3.8802     0.1121  34.606   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 14.36 on 18 degrees of freedom
## Multiple R-squared:  0.9852, Adjusted R-squared:  0.9844 
## F-statistic:  1198 on 1 and 18 DF,  p-value: < 2.2e-16

Koefisien dari $\beta_1$ dan $\beta_2$ dapat kita lihat di bagian coefficients, atau bisa juga dengan menggunakan fungsi coef()

round(coef(lm), 4)

## (Intercept)      income 
##      6.8550      3.8802

Maka kita dapat mengetahui koefisien dari masing-masing $\beta$, sehingga kita dapat menulis estimasi linier regresi kita sebagai berikut:

\[ \begin{align} Y &= \beta_1 + \beta_2 * X \\ Insurance &= 6.855 + 3.8802 * income \end{align} \]

$\beta_2$ merupakan parameter gradien untuk fungsi linier kita. Nilai dari $\beta_2$ menunjukkan bahwa biaya untuk insuransi yang digunakan akan meningkata sebanyak 3.8802 setiap naiknya 1 unit pendapatan rumah tangga.

Kita juga dapat menambahkan garis regresi ke dalam plot kita di awal dengan fungsi abline()

plot(insur$income, insur$insurance,
     ylim = c(0, max(insur$insurance)),
     xlim = c(0, max(insur$income)),
     ylab = "Insurance in $1000",
     xlab = "Income in $1000"
     )
abline(lm)

2.2 Uji Hipotesis

Kita perlu melakukan uji hipotesis untuk mengetahui apakah ada hubungan linier antara variabel bebas dan variabel terikat. Hipotesis kita adalah sebagai berikut:

\[ \begin {align} H_0&: \beta_2 = 0 \\ H_a&: \beta_2 \ne 0 \end{align} \]

Kita akan melakukan tes hipotesis dengan menguji p-value, dimana jika p-value lebih dari tingkat kesalahan, maka kita akan menerima $H_0$ yang berarti tidak ada hubungan linier antara insurance dan income. Sedangkan jika p-value lebih rendah dari tingkat kesalahan, kita akan menolak $H_0$ yang mengartikan bahwa terdapat hubungan linier antara insurance dan income.

Tingkat kesalahan yang digunakan dalam penilitian pada umumnya adalah 5%, sehingga disini saya akan menggunakan tingkat kesalahan 5%. Selanjutnya kita perlu mengecek p-value yang terdapat juga di dalam summary()

as.data.frame(slm$coefficients)

##             Estimate Std. Error    t value     Pr(>|t|)
## (Intercept) 6.854991  7.3834730  0.9284236 3.654714e-01
## income      3.880186  0.1121246 34.6060058 6.391509e-18

Dari hasil di atas, bisa kita lihat bahwa p-value dari income adalah $6.392 \times 10^{-18}$, yang mana lebih rendah dari tingkat kesalahan 0.05 sehingga kita akan menolak $H_0$ yang mengartikan terdapat hubungan linier antara income dan insurance.

2.3 Korelasi

Selanjutnya kita bisa mengecek korelasi antara Insurance dan Income untuk mengetahui seberapa kuatnya hubungan antara Insurance dan Income menggunakan fungsi cor()

round(cor(insur)[1,2],4)

## [1] 0.9926

Sekarang kita akan menginterpretasikan hasil ini menggunakan peringkat korelasi Spearman berikut:

$|r| = 1$ perfect relationship.
$|r|\ge 0.8$ very strong relationship.
$0.6 \le |r| <0.8$ strong relationship.
$0.4 \le |r| <0.6$ moderate relationship.
$0.2 \le |r| <0.4$ weak relationship.
$0 < |r| < 0.2$ very weak relationship.
$|r|=0$ no relationship.

Korelasi antara insurance dan income adalah 0.9926, yang mana mengartikan bahwa ada hubungan yang sangat kuat antara harga insuransi yang digunakan dengan pendapatan rumah tangga.

2.4 Determinasi

Determinasi membantu kita dalam mengetahui seberapa besar variabel bebas menentukan variabel terikat, dalam kasus ini yaitu seberapa besar pendapatan rumah tangga dapat menentukan nilai dari harga insurasnsi yang digunakan. Determinasi bisa kita lihat dari adjusted $R^2$ yang ada di summary().

round(slm$adj.r.squared,4)

## [1] 0.9844

Maka, pendapatan dapat menentukan sebanyak 98.44% nilai dari harga insuransi yang digunakan.

3 Analisis Regresi Linier Berganda

Carilah contoh penerapan analisis regresi linear berganda dalam ilmu ekonometrik! (Persentasikan temuan anda)

Disini saya akan menggunakan data andy dari package PoEdata.

data("andy")
datatable(andy,
          caption = htmltools::tags$caption(
            style = 'caption-side: bottom; text-align: center;',
            htmltools::em('Table Big Andy Sales')),
            extensions = 'FixedColumns',
            option = list(scrollX = TRUE, fixedColumns = TRUE)
          )

Saya akan mulai dengan menjelaskan maksud dari masing-masing kolom. Kolom pertama sales menggambarkan pendapatan per bulan perusahaan Big Andy’s Hamburger dalam $1000, Kolom kedua price menggambarkan indeks harga semua produk yang dijual oleh Big Andy dalam $1000, kolom ketiga advert menggambarkan pengeluaran untuk pengiklanan tiap bulan dalam $1000.

Dari data di atas ini kita akan menganalisa hubungan linier dari variabel terikat sales dengan variabel bebas lainnya.

3.1 Hubungan Linier

Pertama-tama kita perlu mengecek apakah ada hubungan linier antara variabel terikat kita dan variabel bebas. Kita mulai dari membuat model liniernya

lm <- lm(sales~., data=andy)
slm <- summary(lm)
slm

## 
## Call:
## lm(formula = sales ~ ., data = andy)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -13.4825  -3.1434  -0.3456   2.8754  11.3049 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept) 118.9136     6.3516  18.722  < 2e-16 ***
## price        -7.9079     1.0960  -7.215 4.42e-10 ***
## advert        1.8626     0.6832   2.726  0.00804 ** 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 4.886 on 72 degrees of freedom
## Multiple R-squared:  0.4483, Adjusted R-squared:  0.4329 
## F-statistic: 29.25 on 2 and 72 DF,  p-value: 5.041e-10

Setelah membuat model linier kita bisa mengambil koefisiennya dengan rumus coef()

round(coef(lm),4)

## (Intercept)       price      advert 
##    118.9136     -7.9079      1.8626

Maka model linier kita adalah sebagai berikut

\[ \begin{align} Y &= \beta_1 + \beta_2*X_1 +\beta_3*X_2 \\ sales &= 118.9136-7.9079*price+1.8626*advert \end{align} \]

Selanjutnya kita membuat hipotesis kita \[ \begin {align} H_0&: \beta_2, \beta_3 = 0 \\ H_a&: \beta_2, \beta_3\ne 0 \end{align} \]

Untuk menguji hipotesis ini, kita akan menggunakan p-value yang bisa kita lihat dari summary() di bagian coefficients

as.data.frame(slm$coefficients)[,c(1,4)]

##               Estimate     Pr(>|t|)
## (Intercept) 118.913610 2.214293e-29
## price        -7.907854 4.423997e-10
## advert        1.862584 8.038182e-03

Dari p-valuenya kita bisa melihat, bahwa semua p-value <0.05, yang mengartikan terdapat hubungan antara variabel sales dengan price, dan advert

3.2 Efek Variabel

Sekarang kita akan mengecek efek dari tiap variabel bebas(price, advert) terhadap variabel terikat (sales). Kita dapat mengetahui efeknya dengan mengecek koefisien estimasnya, jika koefisiennya negatif, maka memilik efek negatif, sebaliknya jika koefisiennya positif, maka akan memberikan efek positif. Namun di bagian ini kita akan mengecek efeknya menggunakan fungsi-fungsi dari package effects

plot(allEffects(lm))

Dari plot di atas kita bisa melihat bahwa indeks harga produk memiliki hubungan negatif dengan pendapatan per bulan, sehingga setiap kali indeks harga dinaikkan, pendapatan per bulan akan relatif mengalami penurunan. Sebaliknya pengeluaran pengiklanan per bulan memiliki efek positif terhadap penjualan per bulan, sehingga setiap kali menaikkan pengeluaran pengiklanan, maka pendapatan per bulan pun akan relati meningkat

3.3 Korelasi.

Sekarang kita akan menghitung kuat lemahnya hubungan antar variabel terikat dan bebas. Untuk mengetahui seberapa kuatnya hubungan linier, kita dapat mengecek korelasi dari sales dengan price, advert. Kita dapat mengecek korelasi menggunakan fungsi cor()

Correlation <-tidy(cor(andy)[2:3,1])
colnames(Correlation)[1:2] <- c("Regressor", "Correlation")
Correlation

## # A tibble: 2 x 2
##   Regressor Correlation
##   <chr>           <dbl>
## 1 price          -0.626
## 2 advert          0.222

Sekarang kita akan menginterpretasikan hasil ini menggunakan peringkat korelasi Spearman berikut:

$|r| = 1$ perfect relationship.
$|r|\ge 0.8$ very strong relationship.
$0.6 \le |r| <0.8$ strong relationship.
$0.4 \le |r| <0.6$ moderate relationship.
$0.2 \le |r| <0.4$ weak relationship.
$0 < |r| < 0.2$ very weak relationship.
$|r|=0$ no relationship.

Dari kriteria di atas, kita dapat menyimpulkan dari hasil korelasi sebelumnya bahwa pendapatan per bulan memiliki hubungan yang kuat dengan indek harga produk, dan memiliki hubungan yang lemah dengan pengeluaran untuk pengiklanan.

3.4 Determinasi

Determinasi memberitahu kita seberapa banyak variabel bebas dapat mengidentifikasi nilai dari variabel terikat sales. Determinasi disimbolkan dengan $R^2$ yang dapat dicari dengan algoritma berikut:

slm$adj.r.squared

## [1] 0.4329316

Dari $R^2$ di atas, kita dapat menyimpulkan bahwa variabel bebas price, adert dapat menentukan 43.2932 % dari nilai sales.

3.5 Akurasi Estimasi

Untuk mencari estimasi akurasi, kita akan melihat dari nilai Mean Absolute Percentage Error (MAPE). Rumus dari MAPE adalah sebagai berikut:

\[MAPE=\frac{\sum_{i=1}^{n}\text{%} e_i}{n}\]

b0 <- coef(lm)[[1]]
b1 <- coef(lm)[[2]]
x1 <- andy$price
b2 <- coef(lm)[[3]]
x2 <- andy$advert
andy$sales_pred <- (b0 +b1*x1 +b2*x2) %>% round(2)
andy$error <- abs(andy$sales-andy$sales_pred) %>% round(2)
andy$'%error' <- round(andy$error/andy$sales*100,2)
total_percentage_error <- sum(andy$`%error`)
MAPE <-  total_percentage_error/nrow(andy)
paste(round(MAPE,3),"%")

## [1] "4.981 %"

Jadi kita mendapatkan bahwa rataan presentase kesalahannya adalah 4.981%, yang mana mengartikan akurasi dari estimasi adalah sekitar $1-MAPE$

paste("The accuracy is ", round(100-MAPE,3),"%")

## [1] "The accuracy is  95.019 %"

Jadi akurasi price, advert dalam mengestimasi sales adalah sekitar 95.019%

Yang mana menunjukkan tingkat akurasi yang sangat tinggi, sehingga mengestimasi pendapatan per bulan berdasarkan 2 variabel itu sangat disarankan.

Selain itu, akurasi estimasi juga dapat dilihat dari plot perbaning berikut.

df <- data.frame(andy$sales, andy$sales_pred)
ggplot(df[1:75,], aes(x=1:75)) + 
  geom_line(aes(y=andy$sales_pred[1:75], colour="Estimated")) +
  geom_point(aes(x=1:75, y=andy$sales_pred[1:75], colour="Estimated")) +
  geom_line(aes(y=andy$sales[1:75], colour="Actual"))+ 
  geom_point(aes(x=1:75, y=andy$sales[1:75], colour="Actual")) +
  scale_colour_manual("", values = c(Estimated="red", Actual="blue")) +
  labs(x="Observation's Number", y="Wage Value")

Dari plot ini kita dapat melihat banyak titik estimasi yang sama dengan actual walaupun ada beberapa yang berbeda dengan nilai aktual, sehingga model linier kita memiliki tingkat akurasi yang cukup tinggi.

3.6 Multikolinieritas

Dalam menganalisa regresi linier berganda, multikolinieritas merupakan salah satu faktor penting yang harus kita perhatikan. Multikolinieritas adalah kejadian dimana terdapat korelasi yang terlalu tinggin antar variabel bebas.

Untuk mengecek apakah terjadi kasus miltikolinieritas, kita bisa mengecek nilai VIF (Variance Inflation Factor) untuk setiap variabel bebas. Di R, kita dapat mengecek nilai VIF hanya dengan menggunakan fungsi vif() yang mana fungsi ini disediakan oleh package car

vif<-tidy(vif(lm))
colnames(vif)[1:2] <- c("Regressor", "VIF")
vif

## # A tibble: 2 x 2
##   Regressor   VIF
##   <chr>     <dbl>
## 1 price      1.00
## 2 advert     1.00

Dari hasil VIF di atas, kita melihat bahwa tidak ada nilai VIF yang lebih dari 5, yang mana menunjukkan bahwa masalah multikolinieritas tidak terjadi.

Jika anda masih ragu, kita juga dapat mengecek ada tidaknya masalah ini berdasarkan nilai korelasi dari fungsi cor()

cor(andy[,-c(1,4:8)])

##             price     advert
## price  1.00000000 0.02636585
## advert 0.02636585 1.00000000

Bisa kita lihat dari hasil di atas, semua nilai korelasi antar variabel kurang dari 0.8, yang mana mengartikan bahwa masalah multikolinieritas benar-benar tidak terjadi.

4 Pemilihan Model Linier Terbaik

Sehubungan dengan soal No 3, buatlah model regresi linear sederhana yang terbaik dari semua kemungkinan variable (coba terapkan semua kemungkinan model, contohnya, kuardatik, log-log, dll sampai anda menemukan model terbaiknya)

Untuk melihat model terbaik kita akan melihat nilai $R^2$, AIC, dan BIC. Model terbaik adalah model yang memiliki $R^2$ tertinggi, dan BIC serta AIC terendah.

lm1 <- lm(sales~price*advert, data=andy)
lm2<- lm(sales~I(price^2)+I(advert^2)+price*advert, data=andy)
lm3<- lm(sales~I(price^2)+I(advert^2)+I(price^3)+I(advert^3)+price*advert+I(price^2)*advert+price*I(advert^2), data=andy)
lm4<- lm(log(sales)~log(price)+log(advert), data=andy)
lm5<- lm(log(sales)~price+advert, data=andy)
lm6<- lm(sales~log(price)+log(advert), data=andy)
r0 <- as.numeric(glance(lm))
r1 <- as.numeric(glance(lm1))
r2 <- as.numeric(glance(lm2))
r3 <- as.numeric(glance(lm3))
r4 <- as.numeric(glance(lm4))
r5 <- as.numeric(glance(lm5))
r6 <- as.numeric(glance(lm6))
tab <- data.frame(rbind(r0, r1, r2, r3, r4, r5, r6))[,c(1,2,8,9)]
colnames(tab) <- c("Rsq","AdjRsq","AIC","BIC")
row.names(tab) <- c("Biasa", "Interaksi", "Kuadratik", "Polinomial", "Log-Log",  "Log-Lin","Lin-Log")
datatable(tab,
          caption = htmltools::tags$caption(
            style = 'caption-side: bottom; text-align: center;',
            htmltools::em('Comparison of several models')),
            extensions = 'FixedColumns',
            option = list(scrollX = TRUE, fixedColumns = TRUE)
          )

Dari hasil di atas kita melihat adanya perbedaan untuk setiap kriteria. $R^2$ emmberikan model lm3 sebagai model terbaik, AIC dan BIC memberikan model lm4. Karena untuk pemilihan kriteria model AIC dan BIC lebih disarankan dibandingkan $R^2$, maka model terbaik kita adalah model Log-Log dengan model linier sebagai berikut: \[ \begin{align} Log(Y) &= \beta_1 + \beta_2Log(X_1) +\beta_3Log(X_2) \\ Log(sales) &= \beta_1 + \beta_2Log(price)+\beta_3Log(advert) \end{align} \]

coef(lm4)

## (Intercept)  log(price) log(advert) 
##  5.31994869 -0.57493603  0.04544036

Sehingga model lengkap liniernya adalah sebagai berikut: \[ \begin{align} Log(sales) &= 5.31994869 -0.57493603Log(price)+0.04544036Log(advert) \end{align} \]

Econometrics

MID TERM TEST

Julian Salomo

March 24, 2021