library(tidyverse) # data manipulation
## Warning: package 'tidyverse' was built under R version 4.4.3
## Warning: package 'readr' was built under R version 4.4.3
## Warning: package 'dplyr' was built under R version 4.4.3
## Warning: package 'forcats' was built under R version 4.4.3
## Warning: package 'lubridate' was built under R version 4.4.3
library(broom) # clustering algorithms
library(datarium) # clustering algorithms & visualization
## Warning: package 'datarium' was built under R version 4.4.3
library(ggplot2)

Data

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

Menampilkan 3 data teratas dan 3 data terbawah pada data marketing

Model Regresi Sederhana

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

Didapat Model Linier dari sales = 8.439 + 0.04754Youtube

Garis Prediksi Regresi

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 = "red", 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'

Grafik tersebut menunjukkan hubungan antara variabel youtube (kemungkinan besar mengacu pada jumlah pengeluaran iklan atau impresi di YouTube) dengan sales (penjualan). Titik-titik hitam mewakili data aktual, sementara garis biru adalah garis regresi linier yang mencoba memodelkan hubungan tersebut. Secara umum, terlihat tren positif: semakin tinggi nilai youtube, semakin tinggi pula nilai sales. Namun, terdapat variasi atau penyebaran data yang cukup besar di sekitar garis regresi (ditunjukkan dengan garis-garis vertikal merah), menandakan bahwa selain youtube, ada faktor lain yang juga mempengaruhi sales sehingga prediksi penjualan berdasarkan youtube saja tidak sepenuhnya akurat.

Garis Prediksi Residual vs Fitted (Linearitas)

plot(model, 1)

Plot Residuals vs Fitted ini digunakan untuk mengevaluasi apakah model regresi linear sederhana (sales terhadap youtube) sudah sesuai. Secara umum, residual tersebar cukup merata di sekitar garis horizontal nol, menunjukkan bahwa asumsi linearitas dan homoskedastisitas (varian residual konstan) relatif terpenuhi. Tidak terlihat pola melengkung atau pola sistematis lain yang mencolok, meskipun ada sedikit penyebaran lebih lebar di bagian fitted value yang lebih tinggi, yang bisa mengindikasikan adanya sedikit heteroskedastisitas. Selain itu, beberapa titik yang jauh dari garis bisa dianggap sebagai outlier, seperti yang ditandai dengan angka (misal 26, 179, 86). Secara keseluruhan, model ini tampak cukup layak, meski mungkin ada peluang untuk perbaikan lebih lanjut.

Garis Residual Scale-Location (Homoskedastisitas)

plot(model, 3)

Plot Scale-Location ini digunakan untuk memeriksa asumsi homoskedastisitas (kesamaan varians residual) dalam model regresi sales terhadap youtube. Idealnya, titik-titik pada grafik ini harus tersebar secara acak dan mendatar di sekitar garis horizontal. Namun, dari grafik ini terlihat pola kenaikan, di mana nilai √(standardized residuals) cenderung meningkat seiring dengan bertambahnya fitted values. Ini mengindikasikan adanya gejala heteroskedastisitas, yaitu varians residual bertambah saat prediksi penjualan semakin tinggi. Dengan kata lain, model ini cenderung kurang stabil untuk nilai sales yang lebih besar, sehingga mungkin perlu dilakukan transformasi data (seperti log) atau mempertimbangkan model regresi yang lebih kompleks.

Breush-Pagan Test for assurance

Secara Manual (Tanpa Packages)

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 di R

library(lmtest)
## Warning: package 'lmtest' was built under R version 4.4.3
## Loading required package: zoo
## Warning: package 'zoo' was built under R version 4.4.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

Berdasarkan output yang Anda berikan, terlihat bahwa baik perhitungan manual maupun penggunaan fungsi bptest() dari package lmtest menghasilkan nilai statistik Breusch-Pagan (BP) dan nilai p (p-value) yang sangat mirip. Secara manual, nilai BP yang diperoleh adalah 48.03797 dengan p-value sebesar 4.180434e-12, sedangkan menggunakan package, nilai BP adalah 48.038 dengan p-value 4.18e-12. Kesamaan yang signifikan pada kedua hasil ini menunjukkan bahwa implementasi perhitungan secara manual sesuai dengan metode yang dienkapsulasi dalam package lmtest. Dengan p-value yang sangat kecil (jauh di bawah tingkat signifikansi umum seperti 0.05), kita dapat menyimpulkan bahwa terdapat bukti statistik yang kuat untuk menolak hipotesis nol homoskedastisitas. Ini berarti varians dari residual model tidak konstan, atau dengan kata lain, terjadi heteroskedastisitas dalam model regresi yang diuji.

Transformasi Model

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

Dalam konteks ini, setelah mendeteksi adanya heteroskedastisitas pada model awal, transformasi logaritmik pada ‘sales’ kemungkinan dilakukan untuk mengatasi masalah tersebut dan berpotensi meningkatkan kualitas model regresi secara keseluruhan dengan memenuhi asumsi-asumsi regresi linear yang lebih baik. Sehingga didapatkan model linier terbaru yaitu Sales = 2.189 + 0.003Youtube.

Breush-Pagan Test Model Linier Setelah Transformasi

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

Didapatkan p-value sebesar 0.274 yang artinya p-value lebih dari 0.05, sehingga model tidak menunjukkan gejala heteroskedastisitas atau model bersifat homoskedastisitas.

Plot QQ (Normalitas)

plot(log.model, 2)

Grafik Q-Q Residuals mengevaluasi normalitas residual model regresi log(sales) terhadap youtube. Sebagian besar titik mengikuti garis lurus, menandakan residual yang cukup normal. Namun, beberapa titik di kiri bawah menyimpang signifikan, menunjukkan adanya residual ekstrem yang perlu diperhatikan karena dapat mengindikasikan pelanggaran asumsi normalitas.

Uji Kolmogorov-Smirnov

Secara Manual

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

Menggunakan Packages di R

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

Uji Kolmogorov-Smirnov manual dan menggunakan package menghasilkan statistik D yang mirip (0.066 vs 0.071). P-value dari uji menggunakan package adalah 0.2686, yang tidak signifikan. Kesimpulannya, berdasarkan uji KS, tidak ada bukti kuat untuk menolak asumsi normalitas residual pada model.

Uji Kolmogorov-Smirnov Model Asli

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

Diketahui p-value yang didapatkan kecil yaitu 2.2e-16 yang artinya kurang dari 0.05, sehingga disimpulkan bahwa data tidak berdistribusi normal.

Autokorelasi

Pemeriksaan visual

plot(log.model, 1)

Pada grafik ini, titik-titik residual tampak tersebar secara acak di sekitar garis nol, meskipun ada beberapa titik ekstrem. Secara umum, tidak terlihat pola sistematis yang kuat yang mengindikasikan adanya autokorelasi atau non-linearitas yang signifikan, namun keberadaan titik-titik ekstrem perlu dipertimbangkan lebih lanjut.

Autokorelasi (Durbin-Watson)

Secara Manual

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

Menggunakan Packages di 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

Uji Durbin-Watson manual (DW = 1.819) dan menggunakan package (DW = 1.81, p-value = 0.088) memberikan hasil yang serupa. Keduanya mengindikasikan tidak adanya bukti signifikan untuk autokorelasi positif dalam residual model.

Autokorelasi (Durbin-Watson) Model Asli

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

Dari hasil uji autokorelasi tersebut didapatkan p-value 0.3213 dimana lebih besar dari 0.05, sehingga dapat disimpulkan bahwa tidak terdapat autokorelasi pada model asli.

Model Regresi Berganda

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

Persamaan model yang terbentuk adalah sales = 3.527 + 0.046youtube + 0.188facebook − 0.001newspaper

Garis Prediksi Residual vs Fitted (Linearitas)

library(reshape2)
## Warning: package 'reshape2' was built under R version 4.4.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'

Ketidakacakan sebaran residual pada plot-plot diagnostik mengisyaratkan potensi pelanggaran asumsi linearitas. Oleh karena itu, diperlukan pengujian statistik yang lebih formal untuk mengonfirmasi temuan visual ini dan mengevaluasi secara objektif tingkat deviasi dari linearitas.

Garis Residual Scale-Location (Homoskedastisitas)

marketing$fit.r <- model.mlr$residuals
marketing$fit.p <- model.mlr$fitted.values
marketing$sqrt.fit.r <- sqrt(abs(marketing$fit.r))
marketing_long <- marketing |>
 melt(measure.vars = c("youtube", "facebook", "newspaper", "fit.p"))
ggplot(marketing_long, aes(value, sqrt.fit.r, group = variable)) +
 geom_point(shape = 1) + 
 geom_smooth(method = "loess") + 
 geom_hline(yintercept = 0, linetype = "dashed") + 
 facet_wrap(~ variable, scales = "free") + 
 labs(x = "Fitted values / Original Variable", y = "Square Root of 
|Residuals|") +
 theme_minimal() +
 ggtitle("Scale-Location Plot")
## `geom_smooth()` using formula = 'y ~ x'

Berdasarkan hasil plot, untuk variabel Youtube, Facebook, dan fit.p terlihat pola residual yang membentuk lengkungan, mengindikasikan adanya pelanggaran asumsi homoskedastisitas. Sementara itu, pada variabel Variable, penyebaran residual tampak relatif stabil. Dengan demikian, secara umum dapat disimpulkan bahwa model mengalami masalah heteroskedastisitas.

Breush-Pagan Test (Heteroskedastisitas)

Secara Manual

error3 <- residuals(model.mlr)^2
model.error <- lm(error3 ~ youtube + facebook + newspaper, data = marketing)
chi_stat <- dim(marketing)[1]*summary(model.error)$r.square
p_value <- 1 - pchisq(chi_stat, df = 3)
cat("BP =",chi_stat, "; p-value =", p_value)
## BP = 5.132872 ; p-value = 0.1623222

Menggunakan Packages di R

library(lmtest)
bptest(model.mlr)
## 
##  studentized Breusch-Pagan test
## 
## data:  model.mlr
## BP = 5.1329, df = 3, p-value = 0.1623

Uji Breusch-Pagan manual (BP = 5.133, p-value = 0.1623) dan menggunakan package memberikan hasil yang sangat serupa. Dengan p-value yang tidak signifikan, kita tidak memiliki cukup bukti untuk menyimpulkan adanya heteroskedastisitas dalam model regresi berganda.

Plot QQ (Normalitas)

plot(model.mlr, 2)

Grafik Q-Q Residuals untuk model regresi berganda mengevaluasi normalitas residual. Sebagian besar titik terletak dekat dengan garis lurus putus-putus, mengindikasikan bahwa distribusi residual mendekati normal. Namun, terdapat beberapa titik, terutama di ekor kiri bawah (misalnya, titik 131), yang menyimpang dari garis. Penyimpangan ini menunjukkan adanya beberapa residual yang lebih ekstrem dari yang diharapkan di bawah distribusi normal, yang perlu dipertimbangkan dalam evaluasi asumsi model.

Uji Kolmogorov_Smirnov

Secara Manual

error4 <- residuals(model.mlr)
e_cdf <- ecdf(error4)
t_cdf <- pnorm(error4, mean = 0, sd = sd(error4))
D <- max(abs(e_cdf(error4) - t_cdf))
print(D)
## [1] 0.1269655

Menggunakan Packages di R

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

Uji Kolmogorov-Smirnov manual (D = 0.127) dan menggunakan package (D = 0.132, p-value = 0.0019) memberikan statistik D yang serupa. Namun, p-value yang sangat kecil dari hasil package secara signifikan menolak hipotesis nol normalitas residual, menyimpulkan bahwa residual model regresi berganda tidak terdistribusi normal.

Autokorelasi (Durbin-Watson)

Secara Manual

n <- dim(marketing)[1]
r <- cor(error[-1], error[-n])
DW <- 2*(1-r)
DW
## [1] 1.819105
## [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

Menggunakan Packages di R

dwtest(model.mlr)
## 
##  Durbin-Watson test
## 
## data:  model.mlr
## DW = 2.0836, p-value = 0.7236
## alternative hypothesis: true autocorrelation is greater than 0

Uji Durbin-Watson manual (DW = 1.819) dan menggunakan package (DW = 2.084, p-value = 0.724) memberikan hasil yang konsisten. Keduanya menunjukkan tidak adanya bukti signifikan untuk autokorelasi positif dalam residual model regresi berganda.

Multicolinierity

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

Terlihat korelasi positif yang cukup kuat antara ‘youtube’ dan ‘sales’ (0.78), serta antara ‘facebook’ dan ‘sales’ (0.58). Korelasi positif yang lebih lemah juga terlihat antara ‘newspaper’ dan ‘sales’ (0.23). Sementara itu, korelasi antar variabel prediktor relatif rendah, yaitu antara ‘youtube’ dan ‘facebook’ (0.05), ‘youtube’ dan ‘newspaper’ (0.06), serta ‘facebook’ dan ‘newspaper’ (0.35). Korelasi antar prediktor yang rendah ini mengindikasikan bahwa masalah multikolinearitas kemungkinan tidak terlalu signifikan dalam model ini.

Variance Inflation Faktor (VIF)

Secara Manual

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

Menggunakan Packages di 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

Perhitungan VIF manual dan menggunakan package car menghasilkan nilai VIF yang sama untuk youtube (1.005), facebook (1.145), dan newspaper (1.145). Nilai VIF yang rendah ini mengindikasikan tidak adanya masalah multikolinearitas yang signifikan antar prediktor dalam model.

Variabel baru (Simulasi Multikolinearitas)

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

Didapatkan Model linier Sales = 3.529 + 0.045Youtube + 0.264Facebook + 0.148Newspaper - 0.149X4

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

Didapatkan nilai VIF yang tinggi untuk X4 yaitu sebesar 266.0315.

Kesimpulan

Kesimpulan dari analisis ini menegaskan bahwa package-package R seperti lmtest dan car menawarkan efisiensi dan kelengkapan output, termasuk p-value, yang mempermudah interpretasi dibandingkan perhitungan manual, meskipun keduanya menghasilkan kesimpulan substansial yang serupa. Perhitungan manual tetap berharga untuk pemahaman konseptual yang mendalam dan fleksibilitas modifikasi, namun memerlukan lebih banyak upaya dan keahlian teknis. Oleh karena itu, penggunaan package sangat direkomendasikan untuk analisis praktis karena kemudahan dan kelengkapannya.