library(tidyverse)
## Warning: package 'tidyverse' was built under R version 4.3.3
## Warning: package 'tidyr' was built under R version 4.3.3
## Warning: package 'readr' was built under R version 4.3.2
## Warning: package 'purrr' was built under R version 4.3.3
## Warning: package 'dplyr' was built under R version 4.3.3
## Warning: package 'stringr' was built under R version 4.3.3
## Warning: package 'forcats' was built under R version 4.3.3
## Warning: package 'lubridate' was built under R version 4.3.3
## ── Attaching core tidyverse packages ──────────────────────── tidyverse 2.0.0 ──
## ✔ dplyr 1.1.4 ✔ readr 2.1.4
## ✔ 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)
## Warning: package 'broom' was built under R version 4.3.3
library(datarium)
## Warning: package 'datarium' was built under R version 4.3.3
library(ggplot2)
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
Pada bagian ini, beberapa package dimuat untuk keperluan analisis data dan visualisasi, serta data marketing dari package datarium diambil untuk dianalisis. Dataset ini berisi anggaran iklan pada media YouTube, Facebook, dan koran, serta angka penjualan. Fungsi head() dan tail() dipakai untuk melihat sebagian data awal dan akhir guna memastikan strukturnya benar.
model <- lm(sales ~ youtube, data = marketing)
model
##
## Call:
## lm(formula = sales ~ youtube, data = marketing)
##
## Coefficients:
## (Intercept) youtube
## 8.43911 0.04754
Kode ini membentuk model regresi linier sederhana untuk memprediksi sales berdasarkan variabel youtube. Hasilnya menunjukkan hubungan positif, di mana kenaikan biaya iklan YouTube berpengaruh terhadap kenaikan penjualan.
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 = "pink", linewidth = 0.3)+
theme_bw()
## `geom_smooth()` using formula = 'y ~ x'
Di bagian ini dilakukan visualisasi antara data asli dengan garis
regresi serta jarak antara nilai aktual dengan nilai prediksi
(residual). Tujuannya untuk melihat seberapa baik model memetakan data
asli.
plot(model, 1)
Grafik ini digunakan untuk memeriksa apakah hubungan antara variabel X
dan Y bersifat linear. Jika titik-titik menyebar tanpa pola khusus, maka
asumsi linearitas terpenuhi.
plot(model, 3)
Grafik ini membantu melihat apakah residual memiliki varians yang
konstan. Jika titik-titik menyebar merata di sekitar garis horizontal,
maka tidak ada masalah heteroskedastisitas.
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
library(lmtest)
## Warning: package 'lmtest' was built under R version 4.3.3
## Loading required package: zoo
## Warning: package 'zoo' was built under R version 4.3.3
##
## 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
Kode ini melakukan uji Breusch-Pagan untuk menguji apakah varians residual konstan. Uji dilakukan secara manual dan juga menggunakan fungsi bptest() dari package lmtest. Jika p-value kecil, maka terdapat heteroskedastisitas.
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
bptest(log.model)
##
## studentized Breusch-Pagan test
##
## data: log.model
## BP = 1.1959, df = 1, p-value = 0.2742
Karena model awal menunjukkan adanya heteroskedastisitas, maka dilakukan transformasi log pada sales. Setelah itu, uji Breusch-Pagan diulang dan hasilnya menunjukkan perbaikan karena p-value menjadi lebih besar.
plot(log.model, 2)
plot(log.model, 3)
error <- residuals(log.model)
e_cdf <- ecdf(error)
t_cdf <- pnorm(error, mean = 0, sd = sd(error))
D <- max(abs(e_cdf(error) - t_cdf))
print(D)
## [1] 0.0658037
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
Bagian ini digunakan untuk memeriksa apakah residual dari model logaritmik menyebar normal. Digunakan QQ-plot dan uji Kolmogorov–Smirnov baik secara manual maupun fungsi ks.test(). Nilai D kecil dan p-value tinggi menunjukkan residual bersifat normal.
plot(log.model, 1)
n <- dim(marketing)[1]
r <- cor(error[-1], error[-n])
DW <- 2*(1-r)
DW
## [1] 1.819105
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
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
Uji ini digunakan untuk melihat apakah residual saling berkorelasi antar waktu atau tidak. DW dihitung manual dan juga dicek dengan dwtest(). Hasil menunjukkan tidak ada autokorelasi karena nilai DW mendekati 2.
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
Model ini menambahkan dua variabel baru sebagai prediktor untuk sales, yaitu facebook dan newspaper. Ini disebut regresi linier berganda karena melibatkan lebih dari satu variabel independen.
library(reshape2)
## Warning: package 'reshape2' was built under R version 4.3.3
##
## 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'
Visualisasi ini dilakukan untuk melihat pola residual terhadap setiap
variabel input. Jika garis loess mendekati horizontal dan tidak
membentuk pola, maka hubungan dalam model bisa dianggap linear.
library(corrplot)
## Warning: package 'corrplot' was built under R version 4.3.3
## corrplot 0.95 loaded
corrplot.mixed(cor(marketing[,1:4]), upper = "square")
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)
data.frame(Youtube = VIF_y,
Facebook = VIF_f,
Newspaper = VIF_n)
## Youtube Facebook Newspaper
## 1 1.004611 1.144952 1.145187
library(car)
## Warning: package 'car' was built under R version 4.3.3
## Loading required package: carData
## Warning: package 'carData' was built under R version 4.3.3
##
## 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
Kode ini menghitung nilai VIF baik secara manual maupun dengan package car untuk mengecek multikolinearitas. Jika semua nilai VIF di bawah 5, maka tidak ada multikolinearitas yang mengganggu antar variabel input.
set.seed(123)
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.52894 0.04554 0.26424 0.14847 -0.14936
m_x4 <- lm(X4 ~ newspaper + youtube + facebook, data = marketing)
R2_x4 <- summary(m_x4)$r.square
1/(1- R2_x4)
## [1] 266.0315
vif(model.update)
## youtube facebook newspaper X4
## 1.01127 24.32257 195.61739 266.03155
Variabel X4 ditambahkan untuk mensimulasikan kasus multikolinearitas tinggi, karena dibentuk dari kombinasi facebook dan newspaper. Nilai VIF yang sangat besar menunjukkan bahwa X4 sangat colinear dan sebaiknya tidak dimasukkan dalam model.
Analisis ini dilakukan untuk menguji apakah model regresi linier memenuhi asumsi-asumsi klasik (CLRM), seperti linearitas, normalitas, homoskedastisitas, dan tidak adanya autokorelasi, agar hasil prediksi dan inferensi menjadi valid. Dimulai dengan membangun model regresi sederhana antara sales dan youtube, ditemukan adanya masalah heteroskedastisitas berdasarkan uji Breusch-Pagan, sehingga dilakukan transformasi log pada variabel sales untuk memperbaiki distribusi residual. Uji normalitas menunjukkan bahwa residual pada model logaritmik telah menyebar normal, dan tidak terdapat autokorelasi berdasarkan uji Durbin-Watson. Selanjutnya, model diperluas menjadi regresi berganda dengan menambahkan facebook dan newspaper, lalu diuji multikolinearitas melalui korelasi dan VIF, yang hasilnya menunjukkan tidak ada masalah berarti. Terakhir, ditambahkan variabel X4 yang sengaja dibuat colinear untuk menunjukkan bahwa nilai VIF tinggi dapat merusak kestabilan model. Keseluruhan proses ini penting untuk memastikan model regresi yang dibangun benar-benar layak digunakan dalam membuat kesimpulan dan prediksi dari data pemasaran.