Sebelum kita melakukan analisis, mari kita memasukkan package yang diperlukan

library(tidyverse)
## Warning: package 'tidyverse' was built under R version 4.4.3
## Warning: package 'ggplot2' was built under R version 4.4.3
## Warning: package 'tidyr' was built under R version 4.4.3
## Warning: package 'readr' was built under R version 4.4.3
## Warning: package 'purrr' was built under R version 4.4.3
## Warning: package 'dplyr' was built under R version 4.4.3
## Warning: package 'stringr' 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
## ── 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)
## Warning: package 'broom' was built under R version 4.4.3
library(datarium)
## Warning: package 'datarium' was built under R version 4.4.3
library(ggplot2)
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
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
library(corrplot)
## Warning: package 'corrplot' was built under R version 4.4.3
## corrplot 0.95 loaded
library(car)
## Warning: package 'car' was built under R version 4.4.3
## Loading required package: carData
## Warning: package 'carData' was built under R version 4.4.3
## 
## Attaching package: 'car'
## 
## The following object is masked from 'package:dplyr':
## 
##     recode
## 
## The following object is masked from 'package:purrr':
## 
##     some

Model Regresi Linier Sederhana

Marilah kita tampilkan data yang akan digunakan.

data("marketing")
head(marketing, 3)
tail(marketing, 3)

Dari data yang didapatkan, berikut adalah versi modelnya.

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

Maka didapat sebuah model persamaan linier sebagai berikut Sales = 8.439 + 0.047Youtube

Marilah kita melihat garis regressi yang telah diprediksi.

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'

Garis-garis vertikal merah yang menunjukkan selisih antara nilai prediksi dan nilai aktual (residual) semakin panjang seiring berjalannya waktu, yang mengindikasikan adanya heteroskedastisitas, yaitu variansi residual yang tidak konstan. Hal ini merupakan pelanggaran terhadap salah satu asumsi dasar dalam regresi linier dan menunjukkan bahwa model cenderung kurang akurat dalam memprediksi nilai penjualan pada periode dengan volume yang lebih tinggi.

Kita juga dapat melihat visualisasi dari liniearitas data dan juga apakah data yang digunakan homogen.

plot(model, 1) #Linearity

plot(model, 3) #Homogeneity

Dari data yang ditampilkan, menunjukkan bahwa data berkemungkinan heteroskedastisitas. Oleh karena itu kita dapat melakukan uji Breusch-Pagan untuk mengkonfirmasi heteroskedastisitas.

Breusch-Pagan Test (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

Breusch-Pagan Test (Function)

bptest(model)
## 
##  studentized Breusch-Pagan test
## 
## data:  model
## BP = 48.038, df = 1, p-value = 4.18e-12

Dari Breusch-Pagan Test yang dilakukan, dapat disimpulkan bahwa data kita memang heteroskedastisitas. Dikarenakan hasil p-value kurang dari 0.05, maka kita akan melakukan transformasi data dengan syntax berikut.

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

Setelah ditransformasi, kita bisa melakukan pengecekan apakah data masih heteroskedastisitas apa sudah menjadi homogen.

plot(log.model, 3)

Akibat transformasi data, sekarang data kita lebih memenuhi asumsi klasik regresi linear, dimana distribusi residual lebih seragam dan juga berkemungkinan homoskedastisitas. Agar lebih yakin, mari kita melakukan Breusch-Pagan Test lagi.

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

Karena p-valuenya sudah melebihi dari 0.05, maka model yang sekarang tidak menunjukkan adanya heteroskedastisitas, atau asumsi homoskedastisitas sudah terpenuhi.
Selanjutnya kita bisa mengecek normalitas dari data menggunakan syntax berikut.

plot(log.model, 2)

Kita juga dapat mengecek normalitass dengan menggunakan Kolmogrov-Smirnov dengan syntax berikut.

Normality (Manual)

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

Normality (Function)

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

Dengan syntax manual kita tidak dapat melihat p-valuenya, namun melalui function r kita dapat melihat nilai p-value, dari kedua syntax, kita dapat mengambil kesimpulan bahwa residual terdistribusi normal.

Kita juga dapat mengecek pada model pertama apakah sebelumnya sudah terdistribusi normal atau belum.

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

Ternyata pada model pertama, data tidak terdistribusi normal, oleh karena itu mentransformasi model adalah langkah yang tepat.

Selanjutnya kita dapat memeriksa apakah model memiliki autokorelasi dengan menggunakan syntax berikut.

plot(log.model, 1)

Ada sedikit lekukan pada garis merah (smooth line), tapi tidak terlalu kuat. Ini mungkin menunjukkan sedikit non-linearitas, oleh karena itu diperlukan pemeriksaan lebih lanjut, seperti menggunakan tes Durbin-Watson.

Durbin-Watson (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

Durbin-Watson (Function)

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

Dari tes Durbin-Watson yang dilakukan, nilai p-value ditemukan melalui function dari r, dimana nilainya lebih dari 0.05. Maka dapat disimpulkan tidak terdapat autokorelasi.

Model Regresi Linier Berganda

Regresi Linier 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

Model persamaan linier yang didapatkan adalah sales = 3.527 + 0.046youtube + 0.188facebook − 0.001newspaper

Dari model tersebut, kita dapat memeriksa liniearitas.

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'

Pada plot-plot diatas terlihat bahwa Asumsi linearitas belum sepenuhnya terpenuhi maka dapat dilakukan uji test lain untuk memastikan linieritas nya.

Setelah linieritas, kita juga dapat memeriksa apakah model ini mengalami heteroskedasititas atau homoskedasititas.

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()
## `geom_smooth()` using formula = 'y ~ x'

Bisa dilihat pada youtube, facebook, dan fit.p, residual makin besar nilainya di ujung-ujung — indikasi adanya heteroskedastisitas. Agar lebih yakin, kita dapat melakukan Breusch-Pagan Test.

Breusch-Pagan Test (Manual)

error2 <- residuals(model.mlr)^2
model.error <- lm(error2 ~ 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

Breusch-Pagan Test (Function)

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

Karena p-value lebih dari 0.05 maka dapat disimpulkan dapat bersifat homosekdastisitas.

Setelah itu, kita dapat memeriksa normalitas dari model.

plot(model.mlr, 2)

Berdasarkan Q-Q plot, sebagian besar residual model mengikuti distribusi normal dengan baik, terutama di bagian tengah, meskipun ada penyimpangan kecil di bagian ekor kiri dan kanan yang menunjukkan kemungkinan keberadaan outlier atau heavy tails. Agar kita tahu, kita dapat lakukan tes Kolmogrov-Smirnov.

Kolmogrov-Smirnvov (Manual)

error <- residuals(model.mlr)
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.1269655

Kolmogrov-Smirnov (Function)

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

Berdasarkan uji Kolmogorov-Smirnov, diperoleh nilai p-value sebesar 0.001887, yang lebih kecil dari 0.05, sehingga dapat disimpulkan bahwa residual tidak berdistribusi normal secara statistik.

Selanjutnya kita akan memeriksa apa terjadi autokorelasi pada model kita.

Autokorelasi dengan Durbin-Watson (Manual)

n <- dim(marketing)[1]
r <- cor(error[-1], error[-n])
DW <- 2*(1-r)
DW
## [1] 2.094184
dU <- 1.788 #n= 200 dan k= 3
if(DW > dU && DW < 4-dU){
cat("No Autocorrelation, DW = ", DW)
} else {
cat("Autocorrelation exists, DW = ", DW)
}
## No Autocorrelation, DW =  2.094184

Autokorelasi dengan Durbin-Watson (Function)

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

Karena nilai p value lebih dari 0.05 dapat disimpulkan bahwa tidak terjadi autocorrelation pada model ini.

Multikolerasi

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

Berdasarkan heatmap korelasi, dapat disimpulkan bahwa terdapat hubungan positif yang cukup kuat antara pengeluaran iklan di platform YouTube dengan jumlah penjualan (sales) dengan nilai korelasi sebesar 0.78. Ini menunjukkan bahwa semakin besar pengeluaran untuk iklan di YouTube, maka penjualan cenderung meningkat.

Selain itu, iklan di Facebook juga memiliki hubungan positif terhadap penjualan dengan nilai korelasi 0.58, meskipun kekuatannya lebih rendah dibandingkan YouTube. Sementara itu, pengeluaran untuk iklan di surat kabar (newspaper) menunjukkan korelasi yang sangat lemah terhadap penjualan (0.23), sehingga dampaknya terhadap peningkatan penjualan relatif kecil.

Secara keseluruhan, YouTube merupakan media iklan yang paling efektif dalam mendorong penjualan dibandingkan dengan Facebook dan surat kabar dalam dataset ini.

Sekarang kita dapat memeriksa Variance Inflation Factor/VIF, VIF ini dipakai buat mengecek multikolinearitas antar variabel independent (variabel X) dalam regresi linear.

VIF (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)

VIF (Function)

vif(model.mlr)
##   youtube  facebook newspaper 
##  1.004611  1.144952  1.145187

Terlihat bahwa nilai vif semua variabel kurang dari 4, mengindikasikan tidak ada multikorelasi yang terjadi.

Penambahan 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)
summary(model.update)
## 
## Call:
## lm(formula = sales ~ youtube + facebook + newspaper + X4, data = marketing)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -10.2145  -0.9580   0.3256   1.4265   3.5258 
## 
## Coefficients:
##              Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  3.528939   0.371577   9.497  < 2e-16 ***
## youtube      0.045543   0.001389  32.780  < 2e-16 ***
## facebook     0.264243   0.039402   6.706  2.1e-10 ***
## newspaper    0.148471   0.076176   1.949   0.0527 .  
## X4          -0.149359   0.075876  -1.968   0.0504 .  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 2.008 on 195 degrees of freedom
## Multiple R-squared:  0.8992, Adjusted R-squared:  0.8971 
## F-statistic: 434.9 on 4 and 195 DF,  p-value: < 2.2e-16
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

Setelah menambahkan variabel baru (X4), VIF dan R-squared menunjukkan bahwa tidak ada masalah multikolinearitas yang signifikan.

Kesimpulan

Dengan menggunakan function langsung dari r, kita sudah dapat gambaran/analisis yang cukup bagus untuk memeriksa model, namun dengan syntax manual, kita dapat hasil/nilai yang lebih rinci, namun harus dilakukannya pengecekan berulang agar syntax yang digunakan sudah benar atau belum.