Memanggil Packages yang Digunakan

Sebelum melakukan analisis kita harus memanggil packages dibawah ini terlebih dahulu:

library(tidyverse) 
## ── Attaching core tidyverse packages ──────────────────────── tidyverse 2.0.0 ──
## ✔ dplyr     1.1.4     ✔ readr     2.1.5
## ✔ forcats   1.0.0     ✔ stringr   1.5.1
## ✔ ggplot2   3.5.2     ✔ tibble    3.2.1
## ✔ lubridate 1.9.4     ✔ tidyr     1.3.1
## ✔ purrr     1.0.4     
## ── Conflicts ────────────────────────────────────────── tidyverse_conflicts() ──
## ✖ dplyr::filter() masks stats::filter()
## ✖ dplyr::lag()    masks stats::lag()
## ℹ Use the conflicted package (<http://conflicted.r-lib.org/>) to force all conflicts to become errors
library(broom) 
library(datarium) 
library(ggplot2)

Import Data

Kita akan menggunakan set data R pemasaran dari paket R datarium, yang berisi dampak dari tiga media periklanan (youtube, facebook, dan koran) terhadap penjualan. Data adalah anggaran iklan dalam ribuan dolar beserta penjualan. Eksperimen periklanan telah diulang sebanyak 200 kali.

data("marketing")
head(marketing, 3)
##   youtube facebook newspaper sales
## 1  276.12    45.36     83.04 26.52
## 2   53.40    47.16     54.12 12.48
## 3   20.64    55.08     83.16 11.16
tail(marketing, 3)
##     youtube facebook newspaper sales
## 198  212.40    11.16      7.68 15.36
## 199  340.32    50.40     79.44 30.60
## 200  278.52    10.32     10.44 16.08

Regresi Linear Sederhana

Model Regresi Linear Sederhana: Sales = β0 + β1Y outube + ε Menggunakan Package R:

model <- lm(sales ~ youtube, data = marketing)
model
## 
## Call:
## lm(formula = sales ~ youtube, data = marketing)
## 
## Coefficients:
## (Intercept)      youtube  
##     8.43911      0.04754

Cara Manual:

packageVersion("ggplot2")
## [1] '3.5.2'
model.diag.metrics <- augment(model)
ggplot(model.diag.metrics, aes(youtube, sales)) +
geom_point() +
stat_smooth(method = lm, se = FALSE) +
geom_segment(aes(xend = youtube, yend = .fitted), color = "pink2", size = 0.3)+
theme_bw()
## Warning: Using `size` aesthetic for lines was deprecated in ggplot2 3.4.0.
## ℹ Please use `linewidth` instead.
## This warning is displayed once every 8 hours.
## Call `lifecycle::last_lifecycle_warnings()` to see where this warning was
## generated.
## `geom_smooth()` using formula = 'y ~ x'

Linearitas

Untuk mengecek linearitas, kita lihat plot nilai prediksi vs residual.Kalau modelnya linear, residual menyebar acak dan garis merahnya mendatar di sekitar nol. Kalau ada pola melengkung, berarti hubungan antar variabel mungkin tidak linear.

plot(model, 1)

Homogenitas Varians

Plot scale-location digunakan untuk mengecek apakah variansi residual tetap. Kalau titik-titiknya menyebar merata dan garis merah datar, asumsi terpenuhi.Kalau makin melebar atau menyempit, berarti ada masalah heteroskedastisitas.

plot(model, 3)

Breush-Pagan Test

Untuk mengetahui p-Value

Cara Manual:

error2 <- residuals(model)^2
model.error <- lm(error2 ~ marketing$youtube)
chi_stat <- dim(marketing)[1]*summary(model.error)$r.square
p_value <- 1 - pchisq(chi_stat, df = 1)
cat("BP =",chi_stat, "; p-value =", p_value)
## BP = 48.03797 ; p-value = 4.180434e-12

Menggunakan Packages R:

library(lmtest)
## Loading required package: zoo
## 
## Attaching package: 'zoo'
## The following objects are masked from 'package:base':
## 
##     as.Date, as.Date.numeric
bptest(model)
## 
##  studentized Breusch-Pagan test
## 
## data:  model
## BP = 48.038, df = 1, p-value = 4.18e-12
p_value <- 1 - pchisq(chi_stat, df = 1)
cat("BP =",chi_stat, "; p-value =", p_value)
## BP = 48.03797 ; p-value = 4.180434e-12

Transformasi Log Model

Kita akan melakukan transformasi log sedemikian rupa sehingga: log(sales) = β0 + β1youtube + ε

library(lmtest)

log.model <- lm(log(sales) ~ youtube, data = marketing)
log.model
## 
## Call:
## lm(formula = log(sales) ~ youtube, data = marketing)
## 
## Coefficients:
## (Intercept)      youtube  
##    2.189514     0.003156

Memeriksa Homogenitas Varians berdasarkan Visualisasi

plot(log.model, 3)

Memeriksa Homegenitas Varians berdasarkan Breush-Pagan Test

bptest(log.model)
## 
##  studentized Breusch-Pagan test
## 
## data:  log.model
## BP = 1.1959, df = 1, p-value = 0.2742

Uji Normalitas

Plot QQ digunakan untuk mengecek apakah residual berdistribusi normal.
Kalau titik-titiknya mengikuti garis lurus, berarti asumsi normalitas terpenuhi.
Pada contoh ini, titik-titik mengikuti garis, jadi asumsi normalitas bisa diterima.

plot(log.model, 2)

Kita akan menggunakan uji Kolmogorov-Smirnov untuk mengecek normalitas residual.
Di bawah hipotesis nol (\(H_0\)), residual dianggap berdistribusi normal.
Jika hasil uji tidak signifikan, maka \(H_0\) diterima → residual normal.

error <- residuals(log.model)
ecdf_function <- ecdf(error)
t_cdf <- pnorm(error, mean = 0, sd = sd(error))
e_cdf <- ecdf_function(error)
D <- max(abs(e_cdf - t_cdf))
print(D)
## [1] 0.0658037

UJi Normalitas KS-Test

Mari kita lakukan beberapa perbandingan dengan fungsi th’ks.test’ dari R.

Menggunakan Package R (Error):

ks.test(error, "pnorm", 0, sd(error))
## 
##  Asymptotic one-sample Kolmogorov-Smirnov test
## 
## data:  error
## D = 0.070804, p-value = 0.2686
## alternative hypothesis: two-sided

Menggunakan Package R (Model):

ks.test(error, "pnorm", 0, sd(residuals(model)))
## 
##  Asymptotic one-sample Kolmogorov-Smirnov test
## 
## data:  error
## D = 0.45656, p-value < 2.2e-16
## alternative hypothesis: two-sided

Auto Korelasi

Pemeriksaan visual dari uji linearitas bisa digunakan untuk mendeteksi pelanggaran asumsi linearitas. Jika terlihat pola tertentu pada plot residual, berarti asumsi linearitas mungkin tidak terpenuhi.

plot(log.model, 1)

Auto Korelasi Uji Durbin-Watson

Statistik DW digunakan untuk menguji apakah ada autocorrelation pada residual. Nilai DW yang mendekati 2 menunjukkan tidak ada autokorelasi, sedangkan nilai jauh dari 2 menandakan adanya autokorelasi positif atau negatif. DW ≈ 2(1 − ρ)

n <- dim(marketing)[1]
r <- cor(error[-1], error[-n])
DW <- 2*(1-r)
DW
## [1] 1.819105

Nilai dL dan dU untuk dan masing-masing adalah 1,758 dan 1,778. Dan ingat bahwa jika , maka kita gagal menolak hipotesis nol (tidak ada autokorelasi).

dU <- 1.778
if(DW > dU && DW < 4-dU){
cat("No Autocorrelation, DW = ", DW)
} else {
cat("Autocorrelation exists, DW = ", DW)
}
## No Autocorrelation, DW =  1.819105

Menggunakan Packages R:

dwtest(log.model)
## 
##  Durbin-Watson test
## 
## data:  log.model
## DW = 1.81, p-value = 0.08847
## alternative hypothesis: true autocorrelation is greater than 0
dwtest(model)
## 
##  Durbin-Watson test
## 
## data:  model
## DW = 1.9347, p-value = 0.3213
## alternative hypothesis: true autocorrelation is greater than 0

Menambahkan Variabel pada Model

Sekarang, kita akan memasukkan ‘facebook’ dan ‘newspaper’ ke dalam model sehingga: sales = β0 + β1youtube + β2facebook + β3newspaper + ε

model.mlr <- lm(sales ~ youtube + facebook + newspaper, data = marketing)
model.mlr
## 
## Call:
## lm(formula = sales ~ youtube + facebook + newspaper, data = marketing)
## 
## Coefficients:
## (Intercept)      youtube     facebook    newspaper  
##    3.526667     0.045765     0.188530    -0.001037

Linearitas

Kita dapat memeriksa sifat linear hubungan antara DV dan IV dengan beberapa cara. Pertama, kita dapat memplot residual berdasarkan nilai-nilai IV.

library(reshape2)
## 
## Attaching package: 'reshape2'
## The following object is masked from 'package:tidyr':
## 
##     smiths
marketing$fit.r <- model.mlr$residuals
marketing$fit.p <- model.mlr$fitted.values
marketing %>%
melt(measure.vars = c("youtube", "facebook", "newspaper", "fit.p")) %>%
ggplot(aes(value, fit.r, group = variable)) +
geom_point(shape = 1) +
geom_smooth(method = loess) +
geom_hline(yintercept = 0) +
facet_wrap(~ variable, scales = "free")
## `geom_smooth()` using formula = 'y ~ x'

## MultiKorelasi Cara termudah untuk mengukur tingkat multikolinearitas adalah dengan melihat matriks korelasi antara masing-masing variabel.

library(corrplot)
## corrplot 0.95 loaded
corrplot.mixed(cor(marketing[,1:4]), upper = "square")

Varians Inflation Factor (VIF)

Persamaan untuk menghitung faktor inflasi varians untuk prediktor j adalah: VIFj = 1/(1 − Rj^2)

m_youtube <- lm(youtube ~ facebook + newspaper, data = marketing)
m_facebook <- lm(facebook ~ youtube + newspaper, data = marketing)
m_newspaper <- lm(newspaper ~ youtube + facebook, data = marketing)

R2_youtube <- summary(m_youtube)$r.sq
R2_facebook <- summary(m_facebook)$r.sq
R2_newspaper <- summary(m_newspaper)$r.sq
VIF_y <- 1/(1-R2_youtube)
VIF_f <- 1/(1-R2_facebook)
VIF_n <- 1/(1-R2_newspaper)

#Melihat Data Frame
data.frame(Youtube = VIF_y,
Facebook = VIF_f,
Newspaper = VIF_n)
##    Youtube Facebook Newspaper
## 1 1.004611 1.144952  1.145187

Menggunakan Packages R:

library(car)
## Loading required package: carData
## 
## Attaching package: 'car'
## The following object is masked from 'package:dplyr':
## 
##     recode
## The following object is masked from 'package:purrr':
## 
##     some
vif(model.mlr)
##   youtube  facebook newspaper 
##  1.004611  1.144952  1.145187

Penambahan Variabel Kolinear

Kita ingin menambahkan X4 = 0.5 .facebook + newspaper +error

marketing$X4 <- 0.5*marketing$facebook + marketing$newspaper + rnorm(n,0,2)
model.update <- lm(sales ~ youtube + facebook + newspaper + X4, data = marketing)
model.update
## 
## Call:
## lm(formula = sales ~ youtube + facebook + newspaper + X4, data = marketing)
## 
## Coefficients:
## (Intercept)      youtube     facebook    newspaper           X4  
##     3.54958      0.04580      0.22084      0.06425     -0.06556

Memeriksa VIF dari X4 dengan cara manual:

m_x4 <- lm(X4 ~ newspaper + youtube + facebook, data = marketing)
R2_x4 <- summary(m_x4)$r.square
1/(1- R2_x4)
## [1] 207.7772

Memeriksa VIF dari X4 dengan Fungsi R:

vif(model.update)
##    youtube   facebook  newspaper         X4 
##   1.005364  18.537114 153.966986 207.777231