Data marketing dari datarium

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 model Regresi Linier Sederhana dari data marketing menggunakan Package R

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

Regresi linier sederhana (MANUAL)

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

Linieritas

plot(model, 1)

Plot Residuals vs Fitted menunjukkan penyebaran residual terhadap nilai prediksi dari model regresi. Secara umum, residual tersebar cukup acak di sekitar garis nol, yang mengindikasikan bahwa asumsi linearitas dan homoskedastisitas (varian residual konstan) relatif terpenuhi. Namun, terdapat sedikit penyebaran yang melebar pada fitted value yang lebih tinggi, mengisyaratkan kemungkinan adanya sedikit heteroskedastisitas. Tidak ada pola lengkung yang kuat, sehingga model linier tetap layak digunakan.

Visualisasi Variansi Homogenitas

plot(model, 3)

Homogenitas tidak terpenuhi karena sebaran residual semakin melebar seiring kenaikan fitted values. Menunjukan adanya varian residual tidak konstan

Breush-Pagan Test dilakukan untuk mengetahui p-value

#Dengan 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

Dari hasil Breusch-Pagan Test diperoleh nilai statistik BP sebesar 48.03797 dengan p-value sebesar 4.18 × 10⁻¹². Karena p-value jauh lebih kecil dari 0.05, maka dapat disimpulkan bahwa terdapat heteroskedastisitas dalam model regresi ini. Artinya, varians error tidak konstan di seluruh rentang prediksi, sehingga asumsi homoskedastisitas dalam regresi linear sederhana tidak terpenuhi.

Membandingkan hasil Breush-Pagan tes dengan perhitungan manual dan melalui fungsi dalam R dengan Packages R

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

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

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

Hasilnya menunjukkan statistik BP sebesar 1.1959 dengan p-value sebesar 0.2742. Karena p-value lebih besar dari 0.05, maka tidak terdapat cukup bukti untuk menolak hipotesis nol, sehingga dapat disimpulkan bahwa model log sudah memenuhi asumsi homoskedastisitas (varian residual konstan).

Normalitas

plot(log.model, 2)

Homogenitas

plot(log.model, 3)

Normalitas dengan menggunakan Kolmogorov-Smirnov

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

#Normalitas menggunakan fungsi dari 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 dilakukan untuk memeriksa normalitas residual dari model hasil transformasi log. Hasil perhitungan manual menunjukkan nilai statistik D sebesar 0.0658. Dengan menggunakan fungsi ks.test di R, diperoleh p-value sebesar 0.2686. Karena p-value lebih besar dari 0.05, maka tidak ada cukup bukti untuk menolak hipotesis nol, sehingga dapat disimpulkan bahwa residual dari model log berdistribusi normal.

Normalitas dari moedel pertama

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

Autocorrelation

plot(log.model, 1)

Autocorrelation check via Durbin - Watson

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

#Autocorrelation check via Durbin - Watson

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

Let’s do some comparison using ‘dwtest’ in 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

Nilai DW adalah sekitar 1.819, yang berada di antara batas kritis (dU = 1.778) dan (4 - dU = 2.222), sehingga dapat disimpulkan tidak ada autokorelasi. Hasil uji menggunakan fungsi dwtest di R memperkuat kesimpulan ini, dengan DW = 1.81 dan p-value = 0.08847, yang lebih besar dari 0.05, sehingga tidak cukup bukti untuk menyatakan adanya autokorelasi positif. Untuk model model, nilai DW adalah 1.9347 dengan p-value 0.3213, yang juga menunjukkan tidak adanya autokorelasi. Secara keseluruhan, kedua model tidak menunjukkan masalah autokorelasi dalam error-nya.

Let’s add some variables in the model (BERGANDA)

model.mlr <- lm(sales ~ youtube + facebook + newspaper, data = marketing)

Let’s write down the model

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

Linearity

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'

Secara umum, pola residual pada facebook dan newspaper terlihat cukup menyebar acak di sekitar nol, mengindikasikan bahwa asumsi linearitas dan homoskedastisitas cukup terpenuhi. Namun, pada youtube dan fit.p, terlihat adanya pola melengkung (non-linear), yang menunjukkan kemungkinan adanya hubungan non-linear antara variabel-variabel ini dengan responnya. Hal ini mengindikasikan bahwa model mungkin belum sepenuhnya menangkap pola sebenarnya dalam data untuk kedua variabel tersebut.

Homogenitas

plot(model.mlr, 3)

Breush-Pagan Test dilakukan untuk mengetahui p-value

#Dengan cara Manual
error2 <- residuals(model.mlr)^2
model.error <- lm(error2 ~ marketing$youtube + marketing$facebook + marketing$newspaper)
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

Membandingkan hasil Breush-Pangan tes dengan perhitungan manual dan melalui fungsi dalam R dengan packages R

library(lmtest)
bptest(model.mlr)
## 
##  studentized Breusch-Pagan test
## 
## data:  model.mlr
## BP = 5.1329, df = 3, p-value = 0.1623
p_value <- 1 - pchisq(chi_stat, df = 3)
cat("BP =",chi_stat, "; p-value =", p_value)
## BP = 5.132872 ; p-value = 0.1623222

Berdasarkan grafik Scale-Location, distribusi residual terlihat cukup menyebar merata di sekitar garis horizontal, meskipun terdapat sedikit pola naik turun, namun tidak terlalu mencolok. Ini mengindikasikan bahwa asumsi homoskedastisitas (varian residual konstan) sebagian besar terpenuhi. Selain itu, hasil uji Breusch-Pagan menunjukkan nilai statistik BP sebesar 5.132872 dengan p-value sebesar 0.162322. Karena p-value lebih besar dari 0.05, maka tidak ada cukup bukti untuk menolak hipotesis nol, sehingga dapat disimpulkan bahwa tidak terdapat masalah heteroskedastisitas dalam model.

Normalitas dengan menggunakan Kolmogorov-Smirnov

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

Normalitas menggunakan fungsi dari R

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 hasil uji normalitas menggunakan Kolmogorov-Smirnov, diperoleh nilai statistik D sebesar 0.13197 dengan p-value sebesar 0.001887. Karena p-value lebih kecil dari 0.05, maka kita menolak hipotesis nol yang menyatakan bahwa residual berdistribusi normal. Dengan demikian, dapat disimpulkan bahwa residual model tidak mengikuti distribusi normal, yang berarti asumsi normalitas error belum terpenuhi.

Autocorrelation

plot(model.mlr, 1)

Autocorrelation check via Durbin - Watson

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

Autocorrelation check via Durbin - Watson

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

Let’s do some comparison using ‘dwtest’ in 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

Berdasarkan hasil uji autokorelasi menggunakan Durbin-Watson, diperoleh nilai DW sebesar 2.094184. Karena nilai ini berada di antara batas aman (lebih besar dari dU = 1.778 dan kurang dari 4-dU), dapat disimpulkan bahwa tidak terdapat autokorelasi pada residual model. Hasil ini juga dikonfirmasi oleh uji dwtest dari paket R, di mana diperoleh DW = 2.0836 dengan p-value sebesar 0.7236. Karena p-value jauh lebih besar dari 0.05, maka tidak ada cukup bukti untuk menyatakan adanya autokorelasi.

Multicolinierity

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")

Variance Inflation Factor (VIF)

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

VIF

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

Let’s compare with ‘vif’ from ‘R’

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

Hasil dari perhitungan manual VIF menggunakan fungsi didapatkan hasil vif(model.mlr) juga menunjukan nilai yang sama persis. Maka perbandingan dari perhitungan manual dan perhitungan menggunakan fungsi yaitu, menghasilkan nilai VIF yang sama, sehingga tidak ada indikasi multikoinearitas (karna semua VIF < 10). Kesimpulan, model aman dari masalah multikoinearitas.

Let’s add some colinear variable

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

Let’s check the VIF of X4

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

compare with vif from car:

vif(model.update)
##   youtube  facebook newspaper        X4 
##   1.01127  24.32257 195.61739 266.03155

Berdasarkan hasil regresi linear dengan sales sebagai variabel dependen dan youtube, facebook, newspaper, serta X4 sebagai variabel independen, diperoleh model dengan nilai R-squared sebesar 0,8992 dan Adjusted R-squared sebesar 0,8971. Ini menunjukkan bahwa sekitar 89,71% variasi pada sales dapat dijelaskan oleh variasi pada keempat variabel tersebut, yang berarti model ini memiliki tingkat penjelasan yang sangat baik. Koefisien untuk youtube adalah 0,0455 dengan nilai p < 0,001, menunjukkan bahwa pengaruhnya terhadap sales signifikan secara statistik. Begitu juga dengan facebook yang memiliki koefisien 0,2642 dan juga sangat signifikan. Untuk newspaper, koefisien sebesar 0,1485 dengan p-value sekitar 0,0527 menunjukkan bahwa pengaruhnya hanya mendekati signifikan (pada tingkat signifikansi 10%). Sementara itu, X4 memiliki koefisien negatif sebesar -0,1494 dengan p-value 0,0504, juga mendekati signifikan. Secara keseluruhan, model ini sangat baik dalam menjelaskan penjualan, dengan pengaruh yang paling besar dan signifikan berasal dari variabel youtube dan facebook, sedangkan newspaper dan X4 memiliki pengaruh yang lebih lemah dan kurang signifikan. Setelah menambahkan variabel baru (X4), VIF dan R-squared menunjukkan bahwa tidak ada masalah multikolinearitas yang signifikan.