Perbandingan Naive Bayes Classifier dengan Bayesian Logistic Regression dalam menentukan kebenaran berita

Alfidhia

2023-03-30

Library

library(tidyverse)
library(caret)
library(GGally)
library(ggplot2)
library(corrplot)
library(bayesplot)
library(skimr)
library(naivebayes)
library(MCMCpack)

Input Data

Data yang digunakan adalah data “news” dengan peubah tipe_berita sebagai peubah respon. Diberikan data tentang informasi yang diperoleh dari 150 artikel-artikel berita online. Di dalam data tersebut terdapat label yang menunjukkan apakah artikel-artikel berita merupakan berita yang asli atau palsu. Data terdiri dari 4 kolom yaitu : (1) kata_dalam_artikel, (2) tanda_seru_judul_artikel, (3) persen_negatif, dan (4) tipe_berita. Berikut adalah deskripsi dari masing masing kolom

kata_dalam_artikel: banyaknya kata dalam artikel
tanda_seru_judul_artikel: apakah judul artikel mengandung tanda seru?
persen_negatif: persentase kata-kata yang memiliki sentimen negatif
tipe_berita: label artikel berita untuk berita asli atau palsu

news <- read.csv("berita_palsu.csv", stringsAsFactors = T)
str(news)

## 'data.frame':    150 obs. of  4 variables:
##  $ kata_dalam_artikel      : int  219 509 494 268 479 220 184 500 677 485 ...
##  $ tanda_seru_judul_artikel: Factor w/ 2 levels "tidak","ya": 1 1 2 1 1 1 1 2 1 1 ...
##  $ persen_negatif          : num  8.47 4.74 3.33 6.09 2.66 3.02 4.1 4.63 2.18 4.22 ...
##  $ tipe_berita             : Factor w/ 2 levels "asli","palsu": 2 1 2 1 2 1 2 2 2 1 ...

Eksplorasi Data

skim(news)

Data summary
Name	news
Number of rows	150
Number of columns	4
_______________________
Column type frequency:
factor	2
numeric	2
________________________
Group variables	None

Variable type: factor

skim_variable	n_missing	complete_rate	ordered	n_unique	top_counts
tanda_seru_judul_artikel	0	1	FALSE	2	tid: 132, ya: 18
tipe_berita	0	1	FALSE	2	asl: 90, pal: 60

Variable type: numeric

skim_variable	n_missing	complete_rate	mean	sd	p0	p25	p50	p75	p100	hist
kata_dalam_artikel	0	1	533.35	596.44	5	253.00	389.00	534.00	5392.00	▇▁▁▁▁
persen_negatif	0	1	3.13	1.36	0	2.21	3.09	4.04	8.47	▃▇▆▁▁

Dari output di atas, dapat terlihat bahwa terdapat 2 data numeric dan 2 data categoric yang dijadikan class factor pada R. Selain itu, dapat terlihat juga bahwa tidak terdapat missing value pada semua data, dengan masing-masing peubah factor memiliki 2 kategori.

barchart(table(news$tanda_seru_judul_artikel))

barchart(table(news$tipe_berita))

prop.table(table(news$tipe_berita))

## 
##  asli palsu 
##   0.6   0.4

Dari plot di atas dapat terlihat bahwa sebagian besar data di news menjawab merupakan berita asli, perbandingannya yaitu 90:60 atau 60% berita merupakan berita asli. Artinya, data sudah cukup balance, tidak perlu dibalancing.

library(DataExplorer)
plot_scatterplot(data = news,by = "kata_dalam_artikel",
                 geom_point_args = list(color="steelblue"),
                 ggtheme = theme_classic(base_size = 10))

Dari hasil di atas, dapat terlihat bahwa tidak terdapat hubungan yang signifikan antara persen_negatif dengan kata_dalam_artikel. Hal ini tentunya baik untuk pemodelan untuk menghindari multikolinearitas.

plot_correlation(data = news,type = "continuous",cor_args = list(method="spearman"))

Dapat terlihat juga bahwa hubungan kedua peubah numeric sangat kecil.

Splitting Data

set.seed(004)
train_index <- sample(1:nrow(news), 0.8 * nrow(news))
train_data <- news[train_index, ]
test_data <- news[-train_index, ]

Naive Bayes Classifier Modelling

Digunakan package naivebayes untuk kasus ini.

library(naivebayes)

Naive bayes memiliki beberapa hyperparameter yang dapat digunakan. Dalam hal ini, digunakan laplace=0 (default parameter) karena pada data sebelumnya tampak tidak akan terdapat masalah perihal “zero probability problem”. Kemudian, usekernel=FALSE dan usepoisson=FALSE (default parameter) menunjukan bahwa model naive bayes yang digunakan adalah Gaussian-Categorical naive bayes. Artinya, asumsi sebaran yang digunakan adalah sebaran normal untuk setiap peubah prediktor numerik.

Jika digunakan usekernel = TRUE, maka model naive bayes akan mengasumsikan sebaran pada data kontinu sesuai dengan sebaran diperoleh dari pemulusan KDE atau Kernel Density Estimation. Dengan menjadikan argumen usekernel = FALSE, maka sebaran yang digunakan pada data kontinu adalah sebaran normal.
Jika digunakan usepoisson = TRUE, maka model naive bayes akan mengasumsikan sebaran pada data diskrit sesuai dengan sebaran poisson. Dengan menjadikan argumen usepoisson = FALSE, maka sebaran yang digunakan pada data diskrit adalah sebaran normal.

Asumsi ini dibuat karena sebaran normal sering digunakan sebagai asumsi awal dalam model statistik karena sifatnya yang simetris dan mudah diinterpretasikan. Meskipun perlu diperhatikan juga bahwa asumsi ini seharusnya diverifikasi dengan melihat sebaran data untuk memastikan bahwa asumsi tersebut sesuai dengan data yang ada.

Prior bagi peubah respon (y)

Untuk mendapatkan prior \(P(y)\), kita perlu mengetahui berapa banyak berita asli (Y=1) dan berita palsu (Y=0) yang terdapat pada data. Seperti yang tertuliskan pada eksplorasi sebelumnya, atau bisa didapatkan dengan syntax sebagai berikut:

table(train_data$tipe_berita)

## 
##  asli palsu 
##    73    47

Sehingga, prior bisa didapatkan dengan menghitung

\[ P(Y=1)=\frac{73}{73+47}=\frac{73}{120}\approx0.6083 \] dan

\[ P(Y=0)=\frac{47}{73+47}=\frac{73}{120}\approx0.3917 \]

model_nb1

model_nb1 <- naive_bayes(tipe_berita ~ tanda_seru_judul_artikel, data = train_data,
                     laplace = 0, usekernel = FALSE,
                     usepoisson = FALSE)
summary(model_nb1)

## 
## ================================== Naive Bayes ================================== 
##  
## - Call: naive_bayes.formula(formula = tipe_berita ~ tanda_seru_judul_artikel,      data = train_data, laplace = 0, usekernel = FALSE, usepoisson = FALSE) 
## - Laplace: 0 
## - Classes: 2 
## - Samples: 120 
## - Features: 1 
## - Conditional distributions: 
##     - Bernoulli: 1
## - Prior probabilities: 
##     - asli: 0.6083
##     - palsu: 0.3917
## 
## ---------------------------------------------------------------------------------

model_nb2

model_nb2 <- naive_bayes(tipe_berita ~ kata_dalam_artikel, data = train_data,
                     laplace = 0, usekernel = FALSE,
                     usepoisson = FALSE)
model_nb2

## 
## ================================== Naive Bayes ================================== 
##  
##  Call: 
## naive_bayes.formula(formula = tipe_berita ~ kata_dalam_artikel, 
##     data = train_data, laplace = 0, usekernel = FALSE, usepoisson = FALSE)
## 
## --------------------------------------------------------------------------------- 
##  
## Laplace smoothing: 0
## 
## --------------------------------------------------------------------------------- 
##  
##  A priori probabilities: 
## 
##      asli     palsu 
## 0.6083333 0.3916667 
## 
## --------------------------------------------------------------------------------- 
##  
##  Tables: 
## 
## --------------------------------------------------------------------------------- 
##  ::: kata_dalam_artikel (Gaussian) 
## --------------------------------------------------------------------------------- 
##                   
## kata_dalam_artikel     asli    palsu
##               mean 607.8630 446.8298
##               sd   784.1057 301.9390
## 
## ---------------------------------------------------------------------------------

model_nb3

model_nb3 <- naive_bayes(tipe_berita ~ kata_dalam_artikel + persen_negatif, 
                         data = train_data,
                     laplace = 0, usekernel = FALSE,
                     usepoisson = FALSE)
model_nb3

## 
## ================================== Naive Bayes ================================== 
##  
##  Call: 
## naive_bayes.formula(formula = tipe_berita ~ kata_dalam_artikel + 
##     persen_negatif, data = train_data, laplace = 0, usekernel = FALSE, 
##     usepoisson = FALSE)
## 
## --------------------------------------------------------------------------------- 
##  
## Laplace smoothing: 0
## 
## --------------------------------------------------------------------------------- 
##  
##  A priori probabilities: 
## 
##      asli     palsu 
## 0.6083333 0.3916667 
## 
## --------------------------------------------------------------------------------- 
##  
##  Tables: 
## 
## --------------------------------------------------------------------------------- 
##  ::: kata_dalam_artikel (Gaussian) 
## --------------------------------------------------------------------------------- 
##                   
## kata_dalam_artikel     asli    palsu
##               mean 607.8630 446.8298
##               sd   784.1057 301.9390
## 
## --------------------------------------------------------------------------------- 
##  ::: persen_negatif (Gaussian) 
## --------------------------------------------------------------------------------- 
##               
## persen_negatif     asli    palsu
##           mean 2.818630 3.503830
##           sd   1.129036 1.559611
## 
## ---------------------------------------------------------------------------------

model_nb4

model_nb4 <- naive_bayes(tipe_berita ~ kata_dalam_artikel + persen_negatif + tanda_seru_judul_artikel, 
                         data = train_data,
                     laplace = 0, usekernel = FALSE,
                     usepoisson = FALSE)
model_nb4

## 
## ================================== Naive Bayes ================================== 
##  
##  Call: 
## naive_bayes.formula(formula = tipe_berita ~ kata_dalam_artikel + 
##     persen_negatif + tanda_seru_judul_artikel, data = train_data, 
##     laplace = 0, usekernel = FALSE, usepoisson = FALSE)
## 
## --------------------------------------------------------------------------------- 
##  
## Laplace smoothing: 0
## 
## --------------------------------------------------------------------------------- 
##  
##  A priori probabilities: 
## 
##      asli     palsu 
## 0.6083333 0.3916667 
## 
## --------------------------------------------------------------------------------- 
##  
##  Tables: 
## 
## --------------------------------------------------------------------------------- 
##  ::: kata_dalam_artikel (Gaussian) 
## --------------------------------------------------------------------------------- 
##                   
## kata_dalam_artikel     asli    palsu
##               mean 607.8630 446.8298
##               sd   784.1057 301.9390
## 
## --------------------------------------------------------------------------------- 
##  ::: persen_negatif (Gaussian) 
## --------------------------------------------------------------------------------- 
##               
## persen_negatif     asli    palsu
##           mean 2.818630 3.503830
##           sd   1.129036 1.559611
## 
## --------------------------------------------------------------------------------- 
##  ::: tanda_seru_judul_artikel (Bernoulli) 
## --------------------------------------------------------------------------------- 
##                         
## tanda_seru_judul_artikel       asli      palsu
##                    tidak 0.98630137 0.74468085
##                    ya    0.01369863 0.25531915
## 
## ---------------------------------------------------------------------------------

Bayesian Logistic Regression

Package yang digunakan pada kasus ini adalah rstanarm.

library(rstanarm)

Analisis Bayesian memerlukan penentuan sebaran prior f(\(\alpha\)) dan f(\(\beta\)) untuk intersep dan vektor koefisien regresi. Saat menggunakan stan_glm, sebaran ini dapat diatur menggunakan argumen prior_intercept dan prior. Fungsi stan_glm mendukung berbagai sebaran prior, yang dijelaskan dalam dokumentasi rstanarm (help(priors, package = ‘rstanarm’)).

Sebagai contoh, misalkan kita memiliki K prediktor dan percaya — sebelum melihat data — bahwa \(\alpha, \beta_1,...,\beta_K\) memiliki peluang yang sama untukpositif atau negatif, tetapi sangat kecil kemungkinannya untuk jauh dari nol. Kepercayaan ini dapat diwakili oleh sebaran normal dengan rata-rata nol dan ragam yang kecil. Untuk memberikan \(\alpha\) dan setiap \(\beta\) prior ini (dengan standar deviasi 1), dalam panggilan ke stan_glm kita akan menyertakan argumen prior_intercept = normal(0,1) dan prior = normal(0,1).

Jika di sisi lain, kita memiliki kepercayaan a-priori yang lebih sedikit yang menyatakan bahwa parameter akan dekat dengan nol maka kita bisa menggunakan standar deviasi yang lebih besar untuk sebaran normal dan/atau sebaran dengan ekor yang lebih besar daripada normal seperti sebaran Student t. Pada kasus ini, digunakan prior Student t dengan 7 derajat kebebasan dan standar deviasi 2.5, yang merupakan adalah prior default yang wajar ketika koefisien parameter dekat dengan nol tetapi memiliki beberapa peluang menjadi besar.

Fungsi yang digunakan pada tahapan ini adalah student_t yang tersedia pada package rstanarm, argumen yang diisi pada fungsi ini adalah sebagai berikut:

t_prior <- student_t(df = 7, location = 0, scale = 2.5)

yang mana df = 7 berarti derajat bebas dari sebaran-t nya adalah 7, kemudian location = 0 berarti rata-rata dari sebaran-t adalah 0, dan scale = 2.5 berarti standar baku bagi sebaran-t adalah sebesar 2.5.

variable t_prior ini kemudian dimasukkan ke dalam argumen prior = t_prior pada fungsi stan_glm().

model_rl1

t_prior <- student_t(df = 7, location = 0, scale = 2.5)
model_rl1 <- stan_glm(tipe_berita ~ tanda_seru_judul_artikel, data = train_data,
                 family = binomial(link = "logit"), 
                 prior = t_prior, prior_intercept = t_prior,
                 seed = 1401201004, refresh=0, iter = 10000)
summary(model_rl1)

## 
## Model Info:
##  function:     stan_glm
##  family:       binomial [logit]
##  formula:      tipe_berita ~ tanda_seru_judul_artikel
##  algorithm:    sampling
##  sample:       20000 (posterior sample size)
##  priors:       see help('prior_summary')
##  observations: 120
##  predictors:   2
## 
## Estimates:
##                              mean   sd   10%   50%   90%
## (Intercept)                -0.7    0.2 -1.0  -0.7  -0.4 
## tanda_seru_judul_artikelya  3.1    1.0  2.0   3.0   4.4 
## 
## Fit Diagnostics:
##            mean   sd   10%   50%   90%
## mean_PPD 0.4    0.1  0.3   0.4   0.5  
## 
## The mean_ppd is the sample average posterior predictive distribution of the outcome variable (for details see help('summary.stanreg')).
## 
## MCMC diagnostics
##                            mcse Rhat n_eff
## (Intercept)                0.0  1.0  12670
## tanda_seru_judul_artikelya 0.0  1.0   9112
## mean_PPD                   0.0  1.0  15352
## log-posterior              0.0  1.0   7162
## 
## For each parameter, mcse is Monte Carlo standard error, n_eff is a crude measure of effective sample size, and Rhat is the potential scale reduction factor on split chains (at convergence Rhat=1).

model_rl2

model_rl2 <- stan_glm(tipe_berita ~ kata_dalam_artikel, data = train_data,
                 family = binomial(link = "logit"), 
                 prior = t_prior, prior_intercept = t_prior,
                 seed = 1401201004, refresh=0, iter = 10000)
summary(model_rl2)

## 
## Model Info:
##  function:     stan_glm
##  family:       binomial [logit]
##  formula:      tipe_berita ~ kata_dalam_artikel
##  algorithm:    sampling
##  sample:       20000 (posterior sample size)
##  priors:       see help('prior_summary')
##  observations: 120
##  predictors:   2
## 
## Estimates:
##                      mean   sd   10%   50%   90%
## (Intercept)        -0.1    0.3 -0.5  -0.1   0.2 
## kata_dalam_artikel  0.0    0.0  0.0   0.0   0.0 
## 
## Fit Diagnostics:
##            mean   sd   10%   50%   90%
## mean_PPD 0.4    0.1  0.3   0.4   0.5  
## 
## The mean_ppd is the sample average posterior predictive distribution of the outcome variable (for details see help('summary.stanreg')).
## 
## MCMC diagnostics
##                    mcse Rhat n_eff
## (Intercept)        0.0  1.0   7134
## kata_dalam_artikel 0.0  1.0  12212
## mean_PPD           0.0  1.0   5173
## log-posterior      0.0  1.0   4393
## 
## For each parameter, mcse is Monte Carlo standard error, n_eff is a crude measure of effective sample size, and Rhat is the potential scale reduction factor on split chains (at convergence Rhat=1).

model_rl3

model_rl3 <- stan_glm(tipe_berita ~ kata_dalam_artikel + persen_negatif, data = train_data,
                 family = binomial(link = "logit"), 
                 prior = t_prior, prior_intercept = t_prior,
                 QR=T,
                 seed = 1401201004, refresh=0, iter = 10000)
summary(model_rl3)

## 
## Model Info:
##  function:     stan_glm
##  family:       binomial [logit]
##  formula:      tipe_berita ~ kata_dalam_artikel + persen_negatif
##  algorithm:    sampling
##  sample:       20000 (posterior sample size)
##  priors:       see help('prior_summary')
##  observations: 120
##  predictors:   3
## 
## Estimates:
##                      mean   sd   10%   50%   90%
## (Intercept)        -1.4    0.6 -2.1  -1.4  -0.7 
## kata_dalam_artikel  0.0    0.0  0.0   0.0   0.0 
## persen_negatif      0.4    0.2  0.2   0.4   0.6 
## 
## Fit Diagnostics:
##            mean   sd   10%   50%   90%
## mean_PPD 0.4    0.1  0.3   0.4   0.5  
## 
## The mean_ppd is the sample average posterior predictive distribution of the outcome variable (for details see help('summary.stanreg')).
## 
## MCMC diagnostics
##                    mcse Rhat n_eff
## (Intercept)        0.0  1.0  18824
## kata_dalam_artikel 0.0  1.0  14231
## persen_negatif     0.0  1.0  18953
## mean_PPD           0.0  1.0  18618
## log-posterior      0.0  1.0   8969
## 
## For each parameter, mcse is Monte Carlo standard error, n_eff is a crude measure of effective sample size, and Rhat is the potential scale reduction factor on split chains (at convergence Rhat=1).

model_rl4

model_rl4 <- stan_glm(tipe_berita ~ kata_dalam_artikel + persen_negatif + tanda_seru_judul_artikel, data = train_data,
                 family = binomial(link = "logit"),
                 prior = t_prior, prior_intercept = t_prior,
                 QR=T,
                 seed = 1401201004, refresh=0, iter = 10000)
summary(model_rl4)

## 
## Model Info:
##  function:     stan_glm
##  family:       binomial [logit]
##  formula:      tipe_berita ~ kata_dalam_artikel + persen_negatif + tanda_seru_judul_artikel
##  algorithm:    sampling
##  sample:       20000 (posterior sample size)
##  priors:       see help('prior_summary')
##  observations: 120
##  predictors:   4
## 
## Estimates:
##                              mean   sd   10%   50%   90%
## (Intercept)                -1.4    0.6 -2.1  -1.4  -0.7 
## kata_dalam_artikel          0.0    0.0  0.0   0.0   0.0 
## persen_negatif              0.3    0.2  0.1   0.3   0.5 
## tanda_seru_judul_artikelya  2.9    1.0  1.7   2.8   4.2 
## 
## Fit Diagnostics:
##            mean   sd   10%   50%   90%
## mean_PPD 0.4    0.1  0.3   0.4   0.5  
## 
## The mean_ppd is the sample average posterior predictive distribution of the outcome variable (for details see help('summary.stanreg')).
## 
## MCMC diagnostics
##                            mcse Rhat n_eff
## (Intercept)                0.0  1.0  22116
## kata_dalam_artikel         0.0  1.0  18728
## persen_negatif             0.0  1.0  22449
## tanda_seru_judul_artikelya 0.0  1.0  17446
## mean_PPD                   0.0  1.0  19733
## log-posterior              0.0  1.0   8974
## 
## For each parameter, mcse is Monte Carlo standard error, n_eff is a crude measure of effective sample size, and Rhat is the potential scale reduction factor on split chains (at convergence Rhat=1).

Sensitivitas dan Spesifisitas menggunakan data test dengan cut-off peluang 0.5, 0.6 dan 0.7

sensitivitas <- c()
spesifisitas <- c()
for(j in c(0.5,0.6,0.7)){
  for(i in 1:4){
    prediksi <- predict(get(paste0("model_nb",i)), test_data, type="prob")
    resp <- ifelse(prediksi[,1] >= j, "asli", "palsu")
    conf <- confusionMatrix(as.factor(resp),test_data$tipe_berita,positive = "asli")
    sensitivitas <- c(sensitivitas,conf$byClass[1])
    spesifisitas <- c(spesifisitas,conf$byClass[2])
  }
  
  for(i in 1:4){
    pred <-  predict(model_rl1, newdata = test_data, type = "response")
    resp <- ifelse(pred >= j, "asli", "palsu")
    # posterior classification accuracy
    caret::confusionMatrix(as.factor(resp), test_data$tipe_berita, positive="asli")
    sensitivitas <- c(sensitivitas,conf$byClass[1])
    spesifisitas <- c(spesifisitas,conf$byClass[2])
  } 
}
perbandingan <- data.frame(Model=
  rep(c("model_nb1","model_nb2","model_nb3","model_nb4",
        "model_rl1","model_rl2","model_rl3","model_rl4"),3),
  Threshold=rep(c(0.5,0.6,0.7),each=8),
  sensitivitas,
  spesifisitas
)
perbandingan

Untuk menentukan mana model yang paling tinggi sensitifitasnya, yaitu:

perbandingan[which(perbandingan$sensitivitas==max(perbandingan$sensitivitas)),]

Untuk menentukan mana model yang paling tinggi spesisivitasnya, yaitu:

perbandingan[which(perbandingan$spesifisitas==1),]

Kita dapat memilih model mana yang lebih baik tergantung apakah fokus kita terhadap sensitivitas atau spesifisitas

Mari kita lihat lagi perbandingan antar modelnya

perbandingan

Dari hasil di atas, dapat terilhat bahwa model yang memberikan hasil paling seimbang antara sensitivias dan spesivisitas adalah model_rl dengan threshold 0.5, dimana sensitivitas dan spesifisitasnya sebesar 76%. Ketika thresholdnya dinaikan menjadi 0.6 atau 0.7, dapat terlihat bahwa nilai sensitivitas dan spesifisitas menjadi sangat jomplang atau tidak seimbang. Artinya, threshold terbaik adalah 0.5.