Linear Regression vs Bayesian Linear Regression

Alfidhia Rahman NJ

2023-03-30

Input Data

Data yang digunakan adalah data “Bike” dengan peubah Suhu sebagai peubah penjelas dan peubah Pelancong sebagai peubah respon. Data ini terdiri dari 500 amatan dengan 10 amatan teratas sebagaimana yang tertera pada tabel di bawah.

bike <- read.csv("bikes - kuliah mandiri.csv")
colnames(bike) <- c("No", "Suhu_X", "Pelancong_Y")
bike

Soal, Spesifikasi model, dan penjelasannya

\[y=\beta_0+\beta_1 x+\varepsilon\]

di mana y adalah peubah respon, 0 adalah intercept atau slope, 1 adalah koefisien regresi dari peubah penjelas (x), dan ɛ adalah error yang menyebar Normal(0,\(\sigma^2\)).

Model di atas digunakan untuk menyelesaikan tiga kasus yang ada. Penjelasan masing-masing kasus sebagai berikut.

  1. Kasus pertama menyebutkan bahwa sebaran prior bagi 0 dan 1 adalah Uniform.

  2. Kasus kedua menyebutkan bahwa sebaran prior bagi 0 dan 1 adalah Normal(b0,\(s0^2\)) dan Normal(\(b1\),\(s1^2\)) dengan masing-masing parameter sebaran priornya diketahui. Parameter sebaran prior ini akan kami asumsikan sama dengan statistik yang didapatkan oleh kasus pertama, artinya output yang dihasilkan pada kasus pertama akan digunakan sebagai parameter pada kasus kedua.

Eksplorasi Data

library(ggplot2)
## Warning: package 'ggplot2' was built under R version 4.2.3
a <- ggplot(bike, aes(x=bike$Suhu))+
  geom_histogram(color="darkblue", fill="lightblue", bins=30)+
  labs(title="Histogram Peubah suhu",x="Suhu", y = "Count")+
  theme(plot.title = element_text(size = 14L, face = "bold", hjust =  0.5), plot.subtitle = element_text(size = 11L, face = "plain", hjust = 0.5)) + geom_vline(aes(xintercept=mean(bike$Suhu), color="Mean"), linetype="longdash", lwd=1) + geom_vline(aes(xintercept=median(bike$Suhu), color="Median"), linetype="longdash", lwd=1)


b <- ggplot(bike, aes(x=bike$Suhu)) + geom_boxplot(fill='orange', color="darkred", stat = ) +
  labs(title="Boxplot Peubah suhu",x="Suhu")+
  theme(plot.title = element_text(size = 14L, face = "bold", hjust = 0.5), plot.subtitle = element_text(size = 11L, face = "plain", hjust = 0.5))


gridExtra::grid.arrange(a,b)
## Warning: Use of `bike$Suhu` is discouraged.
## ℹ Use `Suhu` instead.
## Warning: Use of `bike$Suhu` is discouraged.
## ℹ Use `Suhu` instead.
## Use of `bike$Suhu` is discouraged.
## ℹ Use `Suhu` instead.
## Use of `bike$Suhu` is discouraged.
## ℹ Use `Suhu` instead.

Terlihat bahwa histogram dan boxplot dari peubah penjelas, yaitu Suhu, yang memberikan informasi mengenai pola sebaran dari data. Dari gambar tersebut terlihat bahwa data cenderung menyebar mengikuti sebaran normal, yang digambarkan dari histogram dan boxplot yang cukup simetris dan dapat dilihat pula bahwa tidak terdapat pencilan pada data.

c <- ggplot(bike, aes(x=bike$Pelancong))+
  geom_histogram(color="darkblue", fill="cyan4", bins=30)+
  labs(title="Histogram Peubah pelancong",x="Jumlah Pelancong", y = "Count")+ theme(plot.title = element_text(size = 14L, face = "bold", hjust = 0.5), plot.subtitle = element_text(size = 11L, face ="plain", hjust = 0.5))+ geom_vline(aes(xintercept=mean(x=bike$Pelancong), color="Mean"),linetype="longdash", lwd=1) + geom_vline(aes(xintercept=median(x=bike$Pelancong), color="Median"), linetype="longdash", lwd=1)


d <- ggplot(bike, aes(x=bike$Pelancong)) +
  geom_boxplot(fill='darkgreen', color="orange", stat = )+
  labs(title="Boxplot Peubah Pelancong",x="Jumlah Pelancong")+
  theme(plot.title = element_text(size = 14L, face = "bold", hjust = 0.5), plot.subtitle = element_text(size = 11L, face = "plain", hjust = 0.5))


gridExtra::grid.arrange(c,d)
## Warning: Use of `bike$Pelancong` is discouraged.
## ℹ Use `Pelancong` instead.
## Use of `bike$Pelancong` is discouraged.
## ℹ Use `Pelancong` instead.
## Use of `bike$Pelancong` is discouraged.
## ℹ Use `Pelancong` instead.
## Use of `bike$Pelancong` is discouraged.
## ℹ Use `Pelancong` instead.

Terlihat bahwa histogram dan boxplot dari peubah respon, yaitu Jumlah Pelancong. Dari gambar tersebut terlihat bahwa peubah respon memiliki sebaran data yang cenderung simetris dan tidak terdapat pencilan pada data.

Regresi Linear Frekuentis

lm=lm(Pelancong_Y~Suhu_X, bike)
summary(lm)
## 
## Call:
## lm(formula = Pelancong_Y ~ Suhu_X, data = bike)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -3564.3 -1009.3  -116.4   959.1  3184.6 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept) -2168.78     363.64  -5.964 4.67e-09 ***
## Suhu_X         89.23       5.67  15.739  < 2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 1288 on 498 degrees of freedom
## Multiple R-squared:  0.3322, Adjusted R-squared:  0.3308 
## F-statistic: 247.7 on 1 and 498 DF,  p-value: < 2.2e-16

Model yang terbentuk:

Pelancong = -2168.78 + 89.23 Satuan Suhu

Artinya dari model yang terbentuk dapat disimpulkan bahwa ketika suhu naik 1 satuan, maka dugaan rataan akan naik sebesar 89.23. Didapatkan informasi bahwa nilai RSE sebesar 1288, nilai ini akan digunakan sebagai sigma pada fungsi di linear bayes kasus 1.

Didapat Residual Standar error sebesar 1288. Nilai ini akan kita gunakan untuk sigma pada prior di regresi bayesian kasus 1.

Visualisasi

library(ggplot2)
beta=coef(lm)
ggplot(data = bike, aes(x = Suhu_X, y = Pelancong_Y)) +
  geom_point(color = "blue") +
  geom_abline(intercept = beta[1], slope = beta[2], lwd = 1)

Regresi Linear Bayesian

Kasus 1

Diketahui bahwa: \(\beta_0\sim flat;\ \beta_1\sim flat\)

Sehingga kita tentukan pada syntax di bawah, argumen slope.prior dan intcpt.prior adalah “flat”. Sehingga, argumen lain dapat diabaikan.

set.seed(04)
library(Bolstad)
## 
## Attaching package: 'Bolstad'
## The following objects are masked from 'package:stats':
## 
##     IQR, sd, var
kasus1 <- bayes.lin.reg(bike$Pelancong_Y, bike$Suhu_X, slope.prior="flat", intcpt.prior="flat", sigma=1288, plot.data=TRUE) #Nilai sigma diketahui berdasarkan standar error pada OLS frekuentis
## Known standard deviation:  1290
##             Posterior Mean Posterior Std. Deviation
##             -------------- ------------------------
## Intercept:  3482           57.601                  
## Slope:      89.23          5.6691

Rata-rata posterior dari intercept adalah 3482, yang merepresentasikan expected values dari parameter intercept dalam distribusi posterior. Standar deviasi posterior dari intercept adalah 57.601, yang merepresentasikan sebaran dari distribusi posterior untuk parameter intercept. Rata-rata posterior dari slope adalah 89.23, yang merepresentasikan nilai ekspektasi dari parameter slope dalam distribusi posterior. Standar deviasi posterior dari slope adalah 5.6691 yang merepresentasikan sebaran dari distribusi posterior untuk parameter slope.

Dari plot di atas, dapat terlihat bahwa sebaran posterior bagi 0 dan 1 dapat terbentuk dari sebaran prior seragam. Hal ini membuktikan bahwa sebaran posterior dapat dibentuk dari sebaran prior apapun, meskipun sesederhana sebaran seragam. Kemudian, dapat terlihat juga bahwa sebagian besar data berada di antara selang kepercayaan bagi y. Oleh karena itu, dapat dikatakan bahwa dugaan model ini telah baik.

Nilai Posterior Mean dan Posterior Std. Deviation untuk intersep dan slope akan digunakan untuk menentukan statistik \(b0\) dan \(b1\) di kasus 2.

Kasus 2

Sigma prior untuk regresi bayesian kasus 2 menyebar \(X\sim \exp(\lambda)\) dengan \(\lambda\) sebagai parameter, yaitu \(E(X)=1/\lambda\). Sehingga, akan menghasilkan

set.seed(04)
kasus2 <- bayes.lin.reg(bike$Pelancong_Y, bike$Suhu_X, slope.prior="normal", intcpt.prior="normal", mb0=89.23, sb0=5.6691, ma0=3482, sa0=57.601, sigma=rexp(n = 1, rate = 1/1288), plot.data=TRUE)
## Known standard deviation:  221
##             Posterior Mean Posterior Std. Deviation
##             -------------- ------------------------
## Intercept:  3482           9.742                   
## Slope:      89.23          0.95881

Rata-rata posterior dari intercept adalah 3482, yang merepresentasikan expected values dari parameter intercept dalam distribusi posterior. Standar deviasi posterior dari intercept adalah 27.818, yang merepresentasikan sebaran dari distribusi posterior untuk parameter intercept. Rata-rata posterior dari slope adalah 89.23, yang merepresentasikan nilai ekspektasi dari parameter slope dalam distribusi posterior. Standar deviasi posterior dari slope adalah 2.7378, yang merepresentasikan sebaran dari distribusi posterior untuk parameter slope.

Dari plot di atas, dapat terlihat bahwa sebaran posterior bagi 0 dan 1 dapat terbentuk dari sebaran prior normal yang merupakan hasil dari kasus 1. Terdapat perbedaan pada sebaran prior dan posterior. Hal ini menunjukan bahwa sebaran sebelumnya telah berhasil diperbaharui dengan memasukkannya terhadap regresi linear bayes yang baru. Dari plot juga dapat terlihat bahwa sebaran yang terbentuk pada posterior lebih ramping yang menandakan bahwa ragam posteriornya lebih rendah daripada ragam prior. Hal ini menunjukan bahwa hasil di kasus kedua lebih baik daripada hasil pada kasus pertama. Kemudian, dapat terlihat juga bahwa sebagian besar data berada di antara selang kepercayaan bagi y. Oleh karena itu, dapat dinyatakan bahwa dugaan model ini telah baik.

Kesimpulan

Dari dua kasus yang diberikan, dapat disimpulkan bahwa dalam analisis regresi linear bayesian, sebaran prior pada parameter model mempengaruhi distribusi posterior parameter. Sebaran prior pada parameter model dapat berupa Uniform atau Normal, dengan nilai parameter prior yang diberikan atau dihitung berdasarkan kasus sebelumnya.

Kasus pertama menghasilkan model dengan parameter intercept 3482 dan slope 89.23, dengan simpangan baku posterior masing-masing sebesar 57.601 dan 5.6691. Hal ini menunjukkan bahwa dalam sebaran prior yang seragam, parameter intersep dan slope memiliki nilai yang konsisten, namun memiliki tingkat variasi yang cukup signifikan.

Kasus kedua mengambil nilai parameter prior dari hasil kasus pertama, dan menghasilkan model dengan parameter intercept 3482 dan slope 89.23, dengan simpangan baku posterior masing-masing sebesar 27.818 dan 2.7378. Hal ini menunjukkan bahwa dengan menggunakan nilai prior yang lebih spesifik, distribusi posterior parameter menjadi lebih terfokus dan memiliki tingkat variasi yang lebih rendah.

Secara keseluruhan, analisis regresi linear bayesian dapat memberikan hasil yang lebih akurat dan konsisten dalam memodelkan data, tergantung pada sebaran prior yang digunakan. Sebaran prior yang lebih spesifik dan kompleks dapat memberikan distribusi posterior parameter yang lebih terfokus dan akurat, sehingga dapat meningkatkan prediksi model terhadap data yang diobservasi.