# ============================================================
# ANALISIS MULTIPLE LINEAR REGRESSION (MLR)
# DATA HOTEL BINTANG INDONESIA TAHUN 2025
# ============================================================


# ============================================================
# 0. PERSIAPAN: LOAD LIBRARY DAN DATA
# ============================================================

library(readxl)
## Warning: package 'readxl' was built under R version 4.5.3
library(dplyr)
## 
## Attaching package: 'dplyr'
## The following objects are masked from 'package:stats':
## 
##     filter, lag
## The following objects are masked from 'package:base':
## 
##     intersect, setdiff, setequal, union
library(car)
## Warning: package 'car' was built under R version 4.5.3
## Loading required package: carData
## Warning: package 'carData' was built under R version 4.5.3
## 
## Attaching package: 'car'
## The following object is masked from 'package:dplyr':
## 
##     recode
library(lmtest)
## Warning: package 'lmtest' was built under R version 4.5.3
## Loading required package: zoo
## Warning: package 'zoo' was built under R version 4.5.3
## 
## Attaching package: 'zoo'
## The following objects are masked from 'package:base':
## 
##     as.Date, as.Date.numeric
library(nortest)
library(MASS)
## Warning: package 'MASS' was built under R version 4.5.3
## 
## Attaching package: 'MASS'
## The following object is masked from 'package:dplyr':
## 
##     select
library(olsrr)
## Warning: package 'olsrr' was built under R version 4.5.3
## 
## Attaching package: 'olsrr'
## The following object is masked from 'package:MASS':
## 
##     cement
## The following object is masked from 'package:datasets':
## 
##     rivers
library(ggplot2)
## Warning: package 'ggplot2' was built under R version 4.5.3
library(corrplot)
## Warning: package 'corrplot' was built under R version 4.5.3
## corrplot 0.95 loaded
library(moments)

set.seed(123)


# Membaca data Excel
data <- read_excel("C:/Users/ACER/Downloads/saibah/anmul/dataset_anmul.xlsx")
# Melihat informasi data
cat("\n", strrep("=", 70), "\n")
## 
##  ======================================================================
cat("DATA YANG DIGUNAKAN\n")
## DATA YANG DIGUNAKAN
cat(strrep("=", 70), "\n")
## ======================================================================
print(head(data))
## # A tibble: 6 × 8
##   Provinsi       Y_Hotel X1_Wisnus X2_PDRB_Kapita X3_Kepadatan X4_Densitas_Jalan
##   <chr>            <dbl>     <dbl>          <dbl>        <dbl>             <dbl>
## 1 Aceh              3553  20096350          45770           99             0.420
## 2 Sumatera Utara   13570  56910262          78310          218             0.553
## 3 Sumatera Barat    5895  22302808          59549          140             0.510
## 4 Riau              9614  23926171         176385           76             0.260
## 5 Jambi             3242  10877483          92786           77             0.236
## 6 Sumatera Sela…    9238  26065516          80663          103             0.223
## # ℹ 2 more variables: `panjang jalan(km2)` <dbl>, `luas wilayah` <dbl>
cat("\nUkuran data:\n")
## 
## Ukuran data:
print(dim(data))
## [1] 38  8
cat("\nNama variabel:\n")
## 
## Nama variabel:
print(names(data))
## [1] "Provinsi"           "Y_Hotel"            "X1_Wisnus"         
## [4] "X2_PDRB_Kapita"     "X3_Kepadatan"       "X4_Densitas_Jalan" 
## [7] "panjang jalan(km2)" "luas wilayah"
cat("\nStruktur data:\n")
## 
## Struktur data:
str(data)
## tibble [38 × 8] (S3: tbl_df/tbl/data.frame)
##  $ Provinsi          : chr [1:38] "Aceh" "Sumatera Utara" "Sumatera Barat" "Riau" ...
##  $ Y_Hotel           : num [1:38] 3553 13570 5895 9614 3242 ...
##  $ X1_Wisnus         : num [1:38] 20096350 56910262 22302808 23926171 10877483 ...
##  $ X2_PDRB_Kapita    : num [1:38] 45770 78310 59549 176385 92786 ...
##  $ X3_Kepadatan      : num [1:38] 99 218 140 76 77 103 106 284 91 268 ...
##  $ X4_Densitas_Jalan : num [1:38] 0.42 0.553 0.51 0.26 0.236 ...
##  $ panjang jalan(km2): num [1:38] 23874 40085 21494 23362 11569 ...
##  $ luas wilayah      : num [1:38] 56835 72438 42108 89901 49023 ...
cat("\nJumlah missing value:\n")
## 
## Jumlah missing value:
print(colSums(is.na(data)))
##           Provinsi            Y_Hotel          X1_Wisnus     X2_PDRB_Kapita 
##                  0                  0                  0                  0 
##       X3_Kepadatan  X4_Densitas_Jalan panjang jalan(km2)       luas wilayah 
##                  0                  0                  0                  0
# ============================================================
# 1. MENYIAPKAN DATA UNTUK ANALISIS MLR
# ============================================================

cat("\n", strrep("=", 70), "\n")
## 
##  ======================================================================
cat("1. DATA UNTUK ANALISIS MLR\n")
## 1. DATA UNTUK ANALISIS MLR
cat(strrep("=", 70), "\n")
## ======================================================================
# Memilih 5 variabel yang digunakan
mlr_data <- data %>%
  dplyr::select(
    Y_Hotel,
    X1_Wisnus,
    X2_PDRB_Kapita,
    X3_Kepadatan,
    X4_Densitas_Jalan
)


cat("\nJumlah observasi dan variabel:\n")
## 
## Jumlah observasi dan variabel:
print(dim(mlr_data))
## [1] 38  5
cat("\nRingkasan data:\n")
## 
## Ringkasan data:
print(summary(mlr_data))
##     Y_Hotel        X1_Wisnus         X2_PDRB_Kapita    X3_Kepadatan     
##  Min.   :   82   Min.   :   379854   Min.   : 19111   Min.   :    6.00  
##  1st Qu.: 1550   1st Qu.:  4348807   1st Qu.: 55017   1st Qu.:   38.75  
##  Median : 5033   Median : 13232428   Median : 73666   Median :  100.00  
##  Mean   :12102   Mean   : 31587699   Mean   : 89420   Mean   :  678.68  
##  3rd Qu.:13312   3rd Qu.: 26477858   3rd Qu.: 90772   3rd Qu.:  255.50  
##  Max.   :84500   Max.   :217197652   Max.   :367687   Max.   :16155.00  
##  X4_Densitas_Jalan
##  Min.   :0.04687  
##  1st Qu.:0.16200  
##  Median :0.36261  
##  Mean   :0.68181  
##  3rd Qu.:0.62431  
##  Max.   :9.83327
cat("\nData MLR:\n")
## 
## Data MLR:
print(mlr_data)
## # A tibble: 38 × 5
##    Y_Hotel X1_Wisnus X2_PDRB_Kapita X3_Kepadatan X4_Densitas_Jalan
##      <dbl>     <dbl>          <dbl>        <dbl>             <dbl>
##  1    3553  20096350          45770           99             0.420
##  2   13570  56910262          78310          218             0.553
##  3    5895  22302808          59549          140             0.510
##  4    9614  23926171         176385           76             0.260
##  5    3242  10877483          92786           77             0.236
##  6    9238  26065516          80663          103             0.223
##  7    1192   7211781          52305          106             0.439
##  8    5008  27090780          55009          284             0.449
##  9    4130   4511219          75323           91             0.356
## 10   12539   4294670         172459          268             0.648
## # ℹ 28 more rows
# ============================================================
# 2. STATISTIK DESKRIPTIF DAN EKSPLORASI DATA
# ============================================================

cat("\n", strrep("=", 70), "\n")
## 
##  ======================================================================
cat("2. STATISTIK DESKRIPTIF DAN EKSPLORASI DATA\n")
## 2. STATISTIK DESKRIPTIF DAN EKSPLORASI DATA
cat(strrep("=", 70), "\n")
## ======================================================================
# Statistik deskriptif
cat("\nStatistik Deskriptif:\n")
## 
## Statistik Deskriptif:
print(summary(mlr_data))
##     Y_Hotel        X1_Wisnus         X2_PDRB_Kapita    X3_Kepadatan     
##  Min.   :   82   Min.   :   379854   Min.   : 19111   Min.   :    6.00  
##  1st Qu.: 1550   1st Qu.:  4348807   1st Qu.: 55017   1st Qu.:   38.75  
##  Median : 5033   Median : 13232428   Median : 73666   Median :  100.00  
##  Mean   :12102   Mean   : 31587699   Mean   : 89420   Mean   :  678.68  
##  3rd Qu.:13312   3rd Qu.: 26477858   3rd Qu.: 90772   3rd Qu.:  255.50  
##  Max.   :84500   Max.   :217197652   Max.   :367687   Max.   :16155.00  
##  X4_Densitas_Jalan
##  Min.   :0.04687  
##  1st Qu.:0.16200  
##  Median :0.36261  
##  Mean   :0.68181  
##  3rd Qu.:0.62431  
##  Max.   :9.83327
# Matriks korelasi
cat("\nMatriks Korelasi:\n")
## 
## Matriks Korelasi:
cor_matrix <- cor(mlr_data)

print(round(cor_matrix, 3))
##                   Y_Hotel X1_Wisnus X2_PDRB_Kapita X3_Kepadatan
## Y_Hotel             1.000     0.694          0.209        0.477
## X1_Wisnus           0.694     1.000          0.068        0.325
## X2_PDRB_Kapita      0.209     0.068          1.000        0.683
## X3_Kepadatan        0.477     0.325          0.683        1.000
## X4_Densitas_Jalan   0.528     0.317          0.654        0.989
##                   X4_Densitas_Jalan
## Y_Hotel                       0.528
## X1_Wisnus                     0.317
## X2_PDRB_Kapita                0.654
## X3_Kepadatan                  0.989
## X4_Densitas_Jalan             1.000
# Visualisasi matriks korelasi
corrplot(
  cor_matrix,
  method = "color",
  type = "upper",
  addCoef.col = "black",
  tl.col = "black",
  tl.srt = 45,
  title = "Matriks Korelasi Antar Variabel",
  mar = c(0, 0, 2, 0)
)

# Scatter plot matrix
pairs(
  mlr_data,
  main = "Scatter Plot Matrix",
  pch = 19
)

# ============================================================
# 3. ESTIMASI MODEL MULTIPLE LINEAR REGRESSION
# ============================================================

cat("\n", strrep("=", 70), "\n")
## 
##  ======================================================================
cat("3. ESTIMASI MODEL MULTIPLE LINEAR REGRESSION\n")
## 3. ESTIMASI MODEL MULTIPLE LINEAR REGRESSION
cat(strrep("=", 70), "\n")
## ======================================================================
# Membentuk model MLR
model <- lm(
  Y_Hotel ~
    X1_Wisnus +
    X2_PDRB_Kapita +
    X3_Kepadatan +
    X4_Densitas_Jalan,
  data = mlr_data
)


# Ringkasan model
cat("\nRingkasan Model:\n")
## 
## Ringkasan Model:
print(summary(model))
## 
## Call:
## lm(formula = Y_Hotel ~ X1_Wisnus + X2_PDRB_Kapita + X3_Kepadatan + 
##     X4_Densitas_Jalan, data = mlr_data)
## 
## Residuals:
##    Min     1Q Median     3Q    Max 
## -15702  -5324  -1104   3352  46176 
## 
## Coefficients:
##                     Estimate Std. Error t value Pr(>|t|)    
## (Intercept)       -7.838e+03  4.812e+03  -1.629 0.112854    
## X1_Wisnus          2.296e-04  3.677e-05   6.246 4.69e-07 ***
## X2_PDRB_Kapita     2.337e-02  3.916e-02   0.597 0.554797    
## X3_Kepadatan      -1.870e+01  4.953e+00  -3.776 0.000632 ***
## X4_Densitas_Jalan  3.416e+04  7.882e+03   4.334 0.000129 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 10770 on 33 degrees of freedom
## Multiple R-squared:  0.7153, Adjusted R-squared:  0.6808 
## F-statistic: 20.72 on 4 and 33 DF,  p-value: 1.275e-08
# Koefisien regresi
cat("\nKoefisien Regresi:\n")
## 
## Koefisien Regresi:
coef_df <- data.frame(
  Variabel = names(coef(model)),
  Koefisien = round(coef(model), 4),
  Std_Error = round(summary(model)$coefficients[, 2], 4),
  t_value = round(summary(model)$coefficients[, 3], 4),
  p_value = round(summary(model)$coefficients[, 4], 4)
)

print(coef_df)
##                            Variabel  Koefisien Std_Error t_value p_value
## (Intercept)             (Intercept) -7837.5029 4811.6515 -1.6289  0.1129
## X1_Wisnus                 X1_Wisnus     0.0002    0.0000  6.2455  0.0000
## X2_PDRB_Kapita       X2_PDRB_Kapita     0.0234    0.0392  0.5967  0.5548
## X3_Kepadatan           X3_Kepadatan   -18.7048    4.9530 -3.7764  0.0006
## X4_Densitas_Jalan X4_Densitas_Jalan 34161.3416 7881.7540  4.3342  0.0001
# Interval kepercayaan 95%
cat("\nInterval Kepercayaan 95%:\n")
## 
## Interval Kepercayaan 95%:
ci <- confint(model, level = 0.95)

print(round(ci, 4))
##                         2.5 %     97.5 %
## (Intercept)       -17626.8816  1951.8757
## X1_Wisnus              0.0002     0.0003
## X2_PDRB_Kapita        -0.0563     0.1030
## X3_Kepadatan         -28.7819    -8.6278
## X4_Densitas_Jalan  18125.7926 50196.8907
# ============================================================
# 4. UJI SIGNIFIKANSI SIMULTAN (UJI F)
# ============================================================

cat("\n", strrep("=", 70), "\n")
## 
##  ======================================================================
cat("4. UJI SIGNIFIKANSI SIMULTAN (UJI F)\n")
## 4. UJI SIGNIFIKANSI SIMULTAN (UJI F)
cat(strrep("=", 70), "\n")
## ======================================================================
# Mengambil informasi F-test
f_stat <- summary(model)$fstatistic

f_value <- f_stat[1]
df1 <- f_stat[2]
df2 <- f_stat[3]

f_pvalue <- pf(
  f_value,
  df1,
  df2,
  lower.tail = FALSE
)


cat("\nHipotesis:\n")
## 
## Hipotesis:
cat("H0: β1 = β2 = β3 = β4 = 0\n")
## H0: β1 = β2 = β3 = β4 = 0
cat("H1: Minimal ada satu βj ≠ 0\n")
## H1: Minimal ada satu βj ≠ 0
cat("\nHasil Uji F:\n")
## 
## Hasil Uji F:
cat("F-statistic =", round(f_value, 4), "\n")
## F-statistic = 20.7243
cat("df1 =", df1, "\n")
## df1 = 4
cat("df2 =", df2, "\n")
## df2 = 33
cat("p-value =", format(f_pvalue, scientific = TRUE), "\n")
## p-value = 1.275116e-08
if (f_pvalue < 0.05) {
  cat("Kesimpulan: Tolak H0 → Model signifikan pada α = 5%\n")
} else {
  cat("Kesimpulan: Gagal menolak H0 → Model tidak signifikan\n")
}
## Kesimpulan: Tolak H0 → Model signifikan pada α = 5%
# ============================================================
# 5. UJI SIGNIFIKANSI PARSIAL (UJI T)
# ============================================================

cat("\n", strrep("=", 70), "\n")
## 
##  ======================================================================
cat("5. UJI SIGNIFIKANSI PARSIAL (UJI T)\n")
## 5. UJI SIGNIFIKANSI PARSIAL (UJI T)
cat(strrep("=", 70), "\n")
## ======================================================================
cat("\nHipotesis untuk setiap variabel:\n")
## 
## Hipotesis untuk setiap variabel:
cat("H0: βj = 0\n")
## H0: βj = 0
cat("H1: βj ≠ 0\n")
## H1: βj ≠ 0
cat("\nHasil Uji T:\n")
## 
## Hasil Uji T:
for (i in 2:nrow(coef_df)) {

  var_name <- coef_df$Variabel[i]
  t_val <- coef_df$t_value[i]
  p_val <- coef_df$p_value[i]

  sig <- ifelse(
    p_val < 0.001, "***",
    ifelse(
      p_val < 0.01, "**",
      ifelse(
        p_val < 0.05, "*",
        ifelse(p_val < 0.1, ".", "ns")
      )
    )
  )

  cat(
    sprintf(
      "%-20s : t = %8.4f, p = %8.6f %s\n",
      var_name,
      t_val,
      p_val,
      sig
    )
  )
}
## X1_Wisnus            : t =   6.2455, p = 0.000000 ***
## X2_PDRB_Kapita       : t =   0.5967, p = 0.554800 ns
## X3_Kepadatan         : t =  -3.7764, p = 0.000600 ***
## X4_Densitas_Jalan    : t =   4.3342, p = 0.000100 ***
cat("\nKeterangan:\n")
## 
## Keterangan:
cat("*** p < 0.001\n")
## *** p < 0.001
cat("**  p < 0.01\n")
## **  p < 0.01
cat("*   p < 0.05\n")
## *   p < 0.05
cat(".   p < 0.10\n")
## .   p < 0.10
cat("ns  tidak signifikan\n")
## ns  tidak signifikan
# ============================================================
# 6. KOEFISIEN DETERMINASI
# ============================================================

cat("\n", strrep("=", 70), "\n")
## 
##  ======================================================================
cat("6. KOEFISIEN DETERMINASI\n")
## 6. KOEFISIEN DETERMINASI
cat(strrep("=", 70), "\n")
## ======================================================================
r_squared <- summary(model)$r.squared
adj_r_squared <- summary(model)$adj.r.squared
resid_se <- summary(model)$sigma


cat("\nR-squared =", round(r_squared, 4), "\n")
## 
## R-squared = 0.7153
cat(
  "Adjusted R-squared =",
  round(adj_r_squared, 4),
  "\n"
)
## Adjusted R-squared = 0.6808
cat(
  "Residual Standard Error =",
  round(resid_se, 4),
  "\n"
)
## Residual Standard Error = 10766.03
cat(
  "Model menjelaskan",
  round(r_squared * 100, 2),
  "% variasi jumlah kamar hotel\n"
)
## Model menjelaskan 71.53 % variasi jumlah kamar hotel
# ============================================================
# 7. UJI ASUMSI RESIDUAL
# ============================================================

cat("\n", strrep("=", 70), "\n")
## 
##  ======================================================================
cat("7. UJI ASUMSI RESIDUAL\n")
## 7. UJI ASUMSI RESIDUAL
cat(strrep("=", 70), "\n")
## ======================================================================
# Residual
residuals_model <- residuals(model)

# Nilai prediksi
fitted_values <- fitted(model)

# Standardized residual
standardized_resid <- rstandard(model)


# ------------------------------------------------------------
# 7.1 UJI NORMALITAS RESIDUAL
# ------------------------------------------------------------

cat("\n--- 7.1 UJI NORMALITAS RESIDUAL ---\n")
## 
## --- 7.1 UJI NORMALITAS RESIDUAL ---
# Shapiro-Wilk
shapiro_test <- shapiro.test(residuals_model)

cat("\nShapiro-Wilk Test:\n")
## 
## Shapiro-Wilk Test:
cat("W =", round(shapiro_test$statistic, 4), "\n")
## W = 0.7806
cat("p-value =", round(shapiro_test$p.value, 6), "\n")
## p-value = 4e-06
cat(
  "Kesimpulan:",
  ifelse(
    shapiro_test$p.value > 0.05,
    "Residual berdistribusi normal",
    "Residual TIDAK berdistribusi normal"
  ),
  "\n"
)
## Kesimpulan: Residual TIDAK berdistribusi normal
# Kolmogorov-Smirnov
ks_test <- ks.test(
  residuals_model,
  "pnorm",
  mean(residuals_model),
  sd(residuals_model)
)

cat("\nKolmogorov-Smirnov Test:\n")
## 
## Kolmogorov-Smirnov Test:
cat("D =", round(ks_test$statistic, 4), "\n")
## D = 0.1881
cat("p-value =", round(ks_test$p.value, 6), "\n")
## p-value = 0.118959
# Anderson-Darling
ad_test <- ad.test(residuals_model)

cat("\nAnderson-Darling Test:\n")
## 
## Anderson-Darling Test:
cat("A =", round(ad_test$statistic, 4), "\n")
## A = 1.7705
cat("p-value =", round(ad_test$p.value, 6), "\n")
## p-value = 0.000127
# Jarque-Bera
jb_test <- jarque.test(residuals_model)

cat("\nJarque-Bera Test:\n")
## 
## Jarque-Bera Test:
cat("JB =", round(jb_test$statistic, 4), "\n")
## JB = 182.8316
cat("p-value =", format(jb_test$p.value, scientific = TRUE), "\n")
## p-value = 0e+00
# ------------------------------------------------------------
# 7.2 UJI HETEROSKEDASTISITAS
# ------------------------------------------------------------

cat("\n--- 7.2 UJI HETEROSKEDASTISITAS ---\n")
## 
## --- 7.2 UJI HETEROSKEDASTISITAS ---
# Breusch-Pagan
bp_test <- bptest(model)

cat("\nBreusch-Pagan Test:\n")
## 
## Breusch-Pagan Test:
cat("BP =", round(bp_test$statistic, 4), "\n")
## BP = 14.6726
cat("p-value =", round(bp_test$p.value, 6), "\n")
## p-value = 0.005431
# Glejser sesuai MLR Hands On
glejser_test <- bptest(
  model,
  ~ fitted(model)
)

cat("\nGlejser Test:\n")
## 
## Glejser Test:
cat("BP =", round(glejser_test$statistic, 4), "\n")
## BP = 5.7508
cat("p-value =", round(glejser_test$p.value, 6), "\n")
## p-value = 0.016481
if (
  bp_test$p.value > 0.05 &
  glejser_test$p.value > 0.05
) {

  cat(
    "Kesimpulan: Tidak terdapat heteroskedastisitas\n"
  )

} else {

  cat(
    "Kesimpulan: Terdapat indikasi heteroskedastisitas\n"
  )
}
## Kesimpulan: Terdapat indikasi heteroskedastisitas
# ------------------------------------------------------------
# 7.3 UJI AUTOKORELASI
# ------------------------------------------------------------

cat("\n--- 7.3 UJI AUTOKORELASI ---\n")
## 
## --- 7.3 UJI AUTOKORELASI ---
# Durbin-Watson
dw_test <- dwtest(model)

cat("\nDurbin-Watson Test:\n")
## 
## Durbin-Watson Test:
cat("DW =", round(dw_test$statistic, 4), "\n")
## DW = 1.4141
cat("p-value =", round(dw_test$p.value, 6), "\n")
## p-value = 0.022876
if (dw_test$p.value > 0.05) {

  cat("Kesimpulan berdasarkan p-value: Tidak terdapat autokorelasi\n")

} else {

  cat("Kesimpulan berdasarkan p-value: Terdapat indikasi autokorelasi\n")
}
## Kesimpulan berdasarkan p-value: Terdapat indikasi autokorelasi
# Breusch-Godfrey
bg_test <- bgtest(
  model,
  order = 1
)

cat("\nBreusch-Godfrey Test:\n")
## 
## Breusch-Godfrey Test:
cat("LM =", round(bg_test$statistic, 4), "\n")
## LM = 3.3211
cat("p-value =", round(bg_test$p.value, 6), "\n")
## p-value = 0.068396
if (bg_test$p.value > 0.05) {

  cat("Kesimpulan: Tidak terdapat autokorelasi\n")

} else {

  cat("Kesimpulan: Terdapat autokorelasi\n")
}
## Kesimpulan: Tidak terdapat autokorelasi
# ------------------------------------------------------------
# 7.4 UJI MULTIKOLINEARITAS
# ------------------------------------------------------------

cat("\n--- 7.4 UJI MULTIKOLINEARITAS ---\n")
## 
## --- 7.4 UJI MULTIKOLINEARITAS ---
# VIF
vif_values <- vif(model)

cat("\nVIF Values:\n")
## 
## VIF Values:
print(round(vif_values, 4))
##         X1_Wisnus    X2_PDRB_Kapita      X3_Kepadatan X4_Densitas_Jalan 
##            1.1854            2.0766           53.2775           48.6452
cat("\nInterpretasi VIF:\n")
## 
## Interpretasi VIF:
for (i in 1:length(vif_values)) {

  vif_val <- vif_values[i]
  var_name <- names(vif_values)[i]

  status <- ifelse(
    vif_val < 5,
    "Tidak ada multikolinearitas",
    ifelse(
      vif_val < 10,
      "Multikolinearitas moderat",
      "Multikolinearitas serius"
    )
  )

  cat(
    sprintf(
      "%-20s : VIF = %8.4f → %s\n",
      var_name,
      vif_val,
      status
    )
  )
}
## X1_Wisnus            : VIF =   1.1854 → Tidak ada multikolinearitas
## X2_PDRB_Kapita       : VIF =   2.0766 → Tidak ada multikolinearitas
## X3_Kepadatan         : VIF =  53.2775 → Multikolinearitas serius
## X4_Densitas_Jalan    : VIF =  48.6452 → Multikolinearitas serius
# Tolerance
tolerance <- 1 / vif_values

cat("\nTolerance Values:\n")
## 
## Tolerance Values:
print(round(tolerance, 4))
##         X1_Wisnus    X2_PDRB_Kapita      X3_Kepadatan X4_Densitas_Jalan 
##            0.8436            0.4816            0.0188            0.0206
# Condition Number
cat(
  "\nCondition Number =",
  round(kappa(model), 4),
  "\n"
)
## 
## Condition Number = 346422501
# ------------------------------------------------------------
# 7.5 UJI LINEARITAS
# ------------------------------------------------------------

cat("\n--- 7.5 UJI LINEARITAS ---\n")
## 
## --- 7.5 UJI LINEARITAS ---
# Ramsey RESET
reset_test <- resettest(
  model,
  power = 2:3,
  type = "fitted"
)

cat("\nRamsey RESET Test:\n")
## 
## Ramsey RESET Test:
cat("F =", round(reset_test$statistic, 4), "\n")
## F = 6.1497
cat("p-value =", round(reset_test$p.value, 6), "\n")
## p-value = 0.005632
if (reset_test$p.value > 0.05) {

  cat(
    "Kesimpulan: Tidak terdapat bukti spesifikasi model linear bermasalah\n"
  )

} else {

  cat(
    "Kesimpulan: Terdapat indikasi spesifikasi model perlu diperiksa\n"
  )
}
## Kesimpulan: Terdapat indikasi spesifikasi model perlu diperiksa
# ============================================================
# 8. DIAGNOSTIK OUTLIER DAN INFLUENTIAL POINTS
# ============================================================

cat("\n", strrep("=", 70), "\n")
## 
##  ======================================================================
cat("8. DIAGNOSTIK OUTLIER DAN INFLUENTIAL POINTS\n")
## 8. DIAGNOSTIK OUTLIER DAN INFLUENTIAL POINTS
cat(strrep("=", 70), "\n")
## ======================================================================
# Jumlah observasi
n <- nrow(mlr_data)


# Cook's Distance
cooks_d <- cooks.distance(model)

cat("\nCook's Distance:\n")
## 
## Cook's Distance:
cat(
  "Nilai maksimum =",
  round(max(cooks_d), 4),
  "\n"
)
## Nilai maksimum = 115.9722
cat(
  "Jumlah observasi dengan Cook's D > 1 =",
  sum(cooks_d > 1),
  "\n"
)
## Jumlah observasi dengan Cook's D > 1 = 2
# Leverage
leverage <- hatvalues(model)

leverage_threshold <- 2 * length(coef(model)) / n

cat("\nLeverage (Hat Values):\n")
## 
## Leverage (Hat Values):
cat(
  "Nilai maksimum =",
  round(max(leverage), 4),
  "\n"
)
## Nilai maksimum = 0.9914
cat(
  "Threshold =",
  round(leverage_threshold, 4),
  "\n"
)
## Threshold = 0.2632
cat(
  "Jumlah observasi di atas threshold =",
  sum(leverage > leverage_threshold),
  "\n"
)
## Jumlah observasi di atas threshold = 4
# Studentized Residual
std_resid <- rstudent(model)

cat("\nStudentized Residuals:\n")
## 
## Studentized Residuals:
cat(
  "Nilai maksimum absolut =",
  round(max(abs(std_resid)), 4),
  "\n"
)
## Nilai maksimum absolut = 14.8075
cat(
  "Jumlah |std_resid| > 3 =",
  sum(abs(std_resid) > 3),
  "\n"
)
## Jumlah |std_resid| > 3 = 1
# DFBETAS
dfbetas_val <- dfbetas(model)

cat("\nDFBETAS:\n")
## 
## DFBETAS:
cat(
  "Jumlah observasi dengan minimal satu |DFBETAS| > 1 =",
  sum(apply(abs(dfbetas_val), 1, function(x) any(x > 1))),
  "\n"
)
## Jumlah observasi dengan minimal satu |DFBETAS| > 1 = 4
# ============================================================
# 9. VISUALISASI DIAGNOSTIK
# ============================================================

cat("\n", strrep("=", 70), "\n")
## 
##  ======================================================================
cat("9. VISUALISASI DIAGNOSTIK\n")
## 9. VISUALISASI DIAGNOSTIK
cat(strrep("=", 70), "\n")
## ======================================================================
par(mfrow = c(2, 2))


# 1. Residual vs Fitted
plot(
  fitted_values,
  residuals_model,
  xlab = "Fitted Values",
  ylab = "Residuals",
  main = "Residual vs Fitted",
  pch = 19
)

abline(
  h = 0,
  lty = 2
)

lines(
  lowess(
    fitted_values,
    residuals_model
  ),
  lwd = 2
)


# 2. Normal Q-Q Plot
qqnorm(
  residuals_model,
  pch = 19,
  main = "Normal Q-Q Plot"
)

qqline(
  residuals_model,
  lwd = 2
)


# 3. Scale-Location
plot(
  fitted_values,
  sqrt(abs(standardized_resid)),
  xlab = "Fitted Values",
  ylab = "√|Standardized Residuals|",
  main = "Scale-Location",
  pch = 19
)

lines(
  lowess(
    fitted_values,
    sqrt(abs(standardized_resid))
  ),
  lwd = 2
)


# 4. Residual vs Leverage
plot(
  leverage,
  standardized_resid,
  xlab = "Leverage",
  ylab = "Standardized Residuals",
  main = "Residual vs Leverage",
  pch = 19
)

abline(
  h = 0,
  lty = 2
)

abline(
  v = leverage_threshold,
  lty = 2
)

par(mfrow = c(1, 1))


# ============================================================
# 10. UJI BOX-COX
# ============================================================

cat("\n", strrep("=", 70), "\n")
## 
##  ======================================================================
cat("10. UJI BOX-COX\n")
## 10. UJI BOX-COX
cat(strrep("=", 70), "\n")
## ======================================================================
bc <- boxcox(
  model,
  lambda = seq(-2, 2, 0.1)
)

lambda_opt <- bc$x[
  which.max(bc$y)
]


cat(
  "\nOptimal Lambda =",
  round(lambda_opt, 4),
  "\n"
)
## 
## Optimal Lambda = 0.2626
cat("\nInterpretasi umum:\n")
## 
## Interpretasi umum:
cat("λ ≈ 1   → tidak perlu transformasi\n")
## λ ≈ 1   → tidak perlu transformasi
cat("λ ≈ 0   → transformasi log\n")
## λ ≈ 0   → transformasi log
cat("λ ≈ 0.5 → transformasi akar\n")
## λ ≈ 0.5 → transformasi akar
# ============================================================
# 11. PEMILIHAN MODEL (STEPWISE)
# ============================================================

cat("\n", strrep("=", 70), "\n")
## 
##  ======================================================================
cat("11. PEMILIHAN MODEL (STEPWISE)\n")
## 11. PEMILIHAN MODEL (STEPWISE)
cat(strrep("=", 70), "\n")
## ======================================================================
step_model <- step(
  model,
  direction = "both",
  trace = 1
)
## Start:  AIC=710.23
## Y_Hotel ~ X1_Wisnus + X2_PDRB_Kapita + X3_Kepadatan + X4_Densitas_Jalan
## 
##                     Df  Sum of Sq        RSS    AIC
## - X2_PDRB_Kapita     1   41265160 3866208038 708.64
## <none>                            3824942878 710.23
## - X3_Kepadatan       1 1653007450 5477950327 721.88
## - X4_Densitas_Jalan  1 2177384368 6002327246 725.36
## - X1_Wisnus          1 4521117994 8346060872 737.88
## 
## Step:  AIC=708.64
## Y_Hotel ~ X1_Wisnus + X3_Kepadatan + X4_Densitas_Jalan
## 
##                     Df  Sum of Sq        RSS    AIC
## <none>                            3866208038 708.64
## + X2_PDRB_Kapita     1   41265160 3824942878 710.23
## - X3_Kepadatan       1 1685558360 5551766399 720.39
## - X4_Densitas_Jalan  1 2152129888 6018337927 723.46
## - X1_Wisnus          1 4573842478 8440050516 736.31
cat("\nModel hasil Stepwise:\n")
## 
## Model hasil Stepwise:
print(summary(step_model))
## 
## Call:
## lm(formula = Y_Hotel ~ X1_Wisnus + X3_Kepadatan + X4_Densitas_Jalan, 
##     data = mlr_data)
## 
## Residuals:
##    Min     1Q Median     3Q    Max 
## -14929  -5786   -384   3500  46547 
## 
## Coefficients:
##                     Estimate Std. Error t value Pr(>|t|)    
## (Intercept)       -5.588e+03  2.962e+03  -1.887 0.067736 .  
## X1_Wisnus          2.245e-04  3.539e-05   6.342 3.11e-07 ***
## X3_Kepadatan      -1.765e+01  4.585e+00  -3.850 0.000497 ***
## X4_Densitas_Jalan  3.312e+04  7.613e+03   4.350 0.000117 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 10660 on 34 degrees of freedom
## Multiple R-squared:  0.7122, Adjusted R-squared:  0.6868 
## F-statistic: 28.04 on 3 and 34 DF,  p-value: 2.588e-09
# ============================================================
# 12. RINGKASAN HASIL UNTUK PPT
# ============================================================

cat("\n", strrep("=", 70), "\n")
## 
##  ======================================================================
cat("12. RINGKASAN HASIL ANALISIS UNTUK PPT\n")
## 12. RINGKASAN HASIL ANALISIS UNTUK PPT
cat(strrep("=", 70), "\n")
## ======================================================================
cat("\n--- MODEL MLR ---\n")
## 
## --- MODEL MLR ---
cat(
  sprintf(
    "Y_Hotel = %.4f + %.6fX1 + %.4fX2 %.4fX3 + %.4fX4\n",
    coef(model)[1],
    coef(model)[2],
    coef(model)[3],
    coef(model)[4],
    coef(model)[5]
  )
)
## Y_Hotel = -7837.5029 + 0.000230X1 + 0.0234X2 -18.7048X3 + 34161.3416X4
cat("\n--- UJI F ---\n")
## 
## --- UJI F ---
cat(
  "F-statistic =",
  round(f_value, 4),
  "\n"
)
## F-statistic = 20.7243
cat(
  "p-value =",
  format(f_pvalue, scientific = TRUE),
  "\n"
)
## p-value = 1.275116e-08
cat("\n--- UJI T ---\n")
## 
## --- UJI T ---
for (i in 2:nrow(coef_df)) {

  cat(
    sprintf(
      "%s : koefisien = %.6f ; p-value = %.6f\n",
      coef_df$Variabel[i],
      coef_df$Koefisien[i],
      coef_df$p_value[i]
    )
  )
}
## X1_Wisnus : koefisien = 0.000200 ; p-value = 0.000000
## X2_PDRB_Kapita : koefisien = 0.023400 ; p-value = 0.554800
## X3_Kepadatan : koefisien = -18.704800 ; p-value = 0.000600
## X4_Densitas_Jalan : koefisien = 34161.341600 ; p-value = 0.000100
cat("\n--- KOEFISIEN DETERMINASI ---\n")
## 
## --- KOEFISIEN DETERMINASI ---
cat(
  "R-squared =",
  round(r_squared, 4),
  "\n"
)
## R-squared = 0.7153
cat(
  "Adjusted R-squared =",
  round(adj_r_squared, 4),
  "\n"
)
## Adjusted R-squared = 0.6808
cat(
  "Persentase variasi yang dijelaskan =",
  round(r_squared * 100, 2),
  "%\n"
)
## Persentase variasi yang dijelaskan = 71.53 %
cat("\n--- UJI ASUMSI ---\n")
## 
## --- UJI ASUMSI ---
cat(
  "Shapiro-Wilk p-value =",
  format(shapiro_test$p.value, scientific = TRUE),
  "\n"
)
## Shapiro-Wilk p-value = 4.30134e-06
cat(
  "Breusch-Pagan p-value =",
  round(bp_test$p.value, 6),
  "\n"
)
## Breusch-Pagan p-value = 0.005431
cat(
  "Durbin-Watson =",
  round(dw_test$statistic, 4),
  "\n"
)
## Durbin-Watson = 1.4141
cat(
  "Breusch-Godfrey p-value =",
  round(bg_test$p.value, 6),
  "\n"
)
## Breusch-Godfrey p-value = 0.068396
cat(
  "VIF maksimum =",
  round(max(vif_values), 4),
  "\n"
)
## VIF maksimum = 53.2775
cat(
  "Ramsey RESET p-value =",
  round(reset_test$p.value, 6),
  "\n"
)
## Ramsey RESET p-value = 0.005632
cat("\n--- DIAGNOSTIK ---\n")
## 
## --- DIAGNOSTIK ---
cat(
  "Cook's D maksimum =",
  round(max(cooks_d), 4),
  "\n"
)
## Cook's D maksimum = 115.9722
cat(
  "Leverage maksimum =",
  round(max(leverage), 4),
  "\n"
)
## Leverage maksimum = 0.9914
cat(
  "Studentized residual maksimum absolut =",
  round(max(abs(std_resid)), 4),
  "\n"
)
## Studentized residual maksimum absolut = 14.8075
cat("\n", strrep("=", 70), "\n")
## 
##  ======================================================================
cat("ANALISIS MLR SELESAI\n")
## ANALISIS MLR SELESAI
cat(strrep("=", 70), "\n")
## ======================================================================

R Markdown

This is an R Markdown document. Markdown is a simple formatting syntax for authoring HTML, PDF, and MS Word documents. For more details on using R Markdown see http://rmarkdown.rstudio.com.

When you click the Knit button a document will be generated that includes both content as well as the output of any embedded R code chunks within the document. You can embed an R code chunk like this:

summary(cars)
##      speed           dist       
##  Min.   : 4.0   Min.   :  2.00  
##  1st Qu.:12.0   1st Qu.: 26.00  
##  Median :15.0   Median : 36.00  
##  Mean   :15.4   Mean   : 42.98  
##  3rd Qu.:19.0   3rd Qu.: 56.00  
##  Max.   :25.0   Max.   :120.00

Including Plots

You can also embed plots, for example:

Note that the echo = FALSE parameter was added to the code chunk to prevent printing of the R code that generated the plot.