Model Regresi Linier Sederhana

Packages yang digunakan:

library(tidyverse) # data manipulation
library(broom)# clustering algorithms
library(datarium)# clustering algorithms & visualization
library(ggplot2)

Menampilkan 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

Membuat model regresi linier sederhana:

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

Model regresi diperoleh: Sales = 8.439 + 0.047Youtube.

Membuat plot residual:

#Plot Residual
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", linewidth = 0.3)+
theme_bw()

Garis merah menunjukkan bahwa semakin besar nilai youtube, semakin besar pula residual, menandakan adanya heteroskedastisitas.

Plot residual vs fitted untuk memeriksa linearitas:

plot(model, 1)

Grafik menunjukkan semakin besar fitted value, semakin besar variasi residual, mengindikasikan heteroskedastisitas.

Plot skala-lokasi untuk memeriksa homoskedastisitas:

plot(model, 3)

Titik-titik residual yang menyebar tidak merata menunjukkan varian residual berubah-ubah (heteroskedastisitas).

Melakukan Breusch-Pagan Test:

#tidak menggunakan 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 in R
library(lmtest)
bptest(model)
## 
##  studentized Breusch-Pagan test
## 
## data:  model
## BP = 48.038, df = 1, p-value = 4.18e-12

Berdasarkan hasil uji kedua code diatas menunjukkan bahwa p-value sangat kecil (jauh di bawah 0.05), sehingga kita dapat menolak H₀ (hipotesis nol) yang artinya terdapat bukti kuat bahwa model mengalami heteroskedastisitas, yaitu varian residual berubah-ubah tergantung pada nilai prediktor.

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

model persamaan linier sebagai berikut Sales = 2.189 + 0.003Youtube

Breush-Pagan Tes untuk model yang telah di transformasi

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

dapat dilihat di p-value telah lebih dari 0.05 artinya model tidak menunjukkan gejala heteroskedastisitas atau model bersifat homoskedastisitas

Plot Normal Q-Q (untuk memeriksa normalitas residual)

plot(log.model, 2)

Dapat dilihat semua titik jatuh kira-kira sepanjang garis, sehingga kita dapat mengasumsikan bahwa data yang telah di transformasi berdistribusi normal.

Kolmogorov-Smirnov Uji

#tidak menggunakan Packages 
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
#menggunakan Packages from 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

Karena nilai p-value lebih dari alfa 0.05 maka dapat disimpulkan data berdistribusi normal

Kolmogorov-Smirnov Uji pada 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

Dapat dilihat Karena nilai p-value kurang dari alfa 0.05 maka dapat disimpulkan data tidak berdistribusi normal

Autocorrelation dengan Durbin - Watson

#tidak menggunakan Function
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 Function From 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

Karena nilia p-value lebih dari 0.05 maka dapat disimpulkan tidak terdapat autocorrelation

Autocorrelation using Durbin - Watson pada model asli

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

Karena nilia p-value lebih dari 0.05 maka dapat disimpulkan tidak terdapat pula autocorrelation pada model asli

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 sales = 3.527 + 0.046youtube + 0.188facebook − 0.001newspaper

Plot residual vs Fitted (untuk memeriksa linearitas)

library(reshape2)
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")

Sebaran data belum sepenuhnya linear, perlu uji tambahan untuk memastikan.

Plot skala-lokasi untuk memeriksa 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")

Breush-Pagan Test untuk memastikan heteroskedastisitas

#tidak menggunakan Packages 
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
#menggunakan Packages in R
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

Plot Normal Q-Q (untuk memeriksa normalitas residual)

plot(model.mlr, 2)

Kolmogorov-Smirnov Uji

#tidak menggunakan Packages 
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
#menggunakan Packages from 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

Autocorrelation dengan Durbin - Watson

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

Karena nilia p value lebih dari 0.05 dapat disimpulkan bahwa tidak terjadi autocorrelation pada model ini  Multicorrelation

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

Menghitung Variance Inflation Factor (VIF)

#tidak menggunakan Packages
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 from R
library(car)
vif(model.mlr)
##   youtube  facebook newspaper 
##  1.004611  1.144952  1.145187

Semua VIF < 4, menunjukkan tidak ada multikolinearitas serius.

Penambahan variabel baru untuk 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 penambahan variabel X4, VIF tetap rendah, menunjukkan multikolinearitas tidak menjadi masalah besar.