# =========================================================
# ANALISIS MULTIPLE LINEAR REGRESSION (MLR)
# JUMLAH KAMAR HOTEL BINTANG DI INDONESIA TAHUN 2025
# =========================================================


# =========================================================
# 1. PERSIAPAN
# =========================================================

rm(list = ls())

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


# =========================================================
# 2. INPUT DATA
# =========================================================

data <- read_excel(file.choose())

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 ...
View(data)
dim(data)
## [1] 38  8
names(data)
## [1] "Provinsi"           "Y_Hotel"            "X1_Wisnus"         
## [4] "X2_PDRB_Kapita"     "X3_Kepadatan"       "X4_Densitas_Jalan" 
## [7] "panjang jalan(km2)" "luas wilayah"
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
summary(data)
##    Provinsi            Y_Hotel        X1_Wisnus         X2_PDRB_Kapita  
##  Length:38          Min.   :   82   Min.   :   379854   Min.   : 19111  
##  Class :character   1st Qu.: 1550   1st Qu.:  4348807   1st Qu.: 55017  
##  Mode  :character   Median : 5033   Median : 13232428   Median : 73666  
##                     Mean   :12102   Mean   : 31587699   Mean   : 89420  
##                     3rd Qu.:13312   3rd Qu.: 26477858   3rd Qu.: 90772  
##                     Max.   :84500   Max.   :217197652   Max.   :367687  
##   X3_Kepadatan      X4_Densitas_Jalan panjang jalan(km2)  luas wilayah     
##  Min.   :    6.00   Min.   :0.04687   Min.   : 4715      Min.   :   661.5  
##  1st Qu.:   38.75   1st Qu.:0.16200   1st Qu.: 6310      1st Qu.: 19754.5  
##  Median :  100.00   Median :0.36261   Median :10044      Median : 43719.1  
##  Mean   :  678.68   Mean   :0.68181   Mean   :14161      Mean   : 49793.0  
##  3rd Qu.:  255.50   3rd Qu.:0.62431   3rd Qu.:18990      3rd Qu.: 61472.5  
##  Max.   :16155.00   Max.   :9.83327   Max.   :42976      Max.   :153443.9
# =========================================================
# 3. MEMILIH VARIABEL YANG DIGUNAKAN
# =========================================================

mlr_data <- data %>%
  dplyr::select(
    Y_Hotel,
    X1_Wisnus,
    X2_PDRB_Kapita,
    X3_Kepadatan,
    X4_Densitas_Jalan
  )

str(mlr_data)
## tibble [38 × 5] (S3: tbl_df/tbl/data.frame)
##  $ 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 ...
View(mlr_data)
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
# =========================================================
# 4. STATISTIK DESKRIPTIF
# =========================================================

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("\nMEAN:\n")
## 
## MEAN:
print(sapply(mlr_data, mean))
##           Y_Hotel         X1_Wisnus    X2_PDRB_Kapita      X3_Kepadatan 
##      1.210234e+04      3.158770e+07      8.941963e+04      6.786842e+02 
## X4_Densitas_Jalan 
##      6.818084e-01
cat("\nSTANDAR DEVIASI:\n")
## 
## STANDAR DEVIASI:
print(sapply(mlr_data, sd))
##           Y_Hotel         X1_Wisnus    X2_PDRB_Kapita      X3_Kepadatan 
##      1.905422e+04      5.241116e+07      6.512340e+04      2.608284e+03 
## X4_Densitas_Jalan 
##      1.566217e+00
cat("\nVARIANS:\n")
## 
## VARIANS:
print(sapply(mlr_data, var))
##           Y_Hotel         X1_Wisnus    X2_PDRB_Kapita      X3_Kepadatan 
##      3.630632e+08      2.746930e+15      4.241058e+09      6.803148e+06 
## X4_Densitas_Jalan 
##      2.453036e+00
# =========================================================
# 5. KORELASI ANTARVARIABEL
# =========================================================

cor_matrix <- cor(mlr_data)

cat("\nMATRIKS KORELASI:\n")
## 
## MATRIKS KORELASI:
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
corrplot(
  cor_matrix,
  method = "number",
  type = "upper",
  tl.col = "black"
)

# =========================================================
# 6. MODEL MULTIPLE LINEAR REGRESSION
# =========================================================

model <- lm(
  Y_Hotel ~ X1_Wisnus +
    X2_PDRB_Kapita +
    X3_Kepadatan +
    X4_Densitas_Jalan,
  data = mlr_data
)

cat("\nHASIL MODEL MLR:\n")
## 
## HASIL MODEL MLR:
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
# =========================================================
# 7. KOEFISIEN REGRESI
# =========================================================

cat("\nKOEFISIEN REGRESI:\n")
## 
## KOEFISIEN REGRESI:
print(coef(summary(model)))
##                        Estimate   Std. Error    t value     Pr(>|t|)
## (Intercept)       -7.837503e+03 4.811652e+03 -1.6288592 1.128536e-01
## X1_Wisnus          2.296282e-04 3.676696e-05  6.2455029 4.692649e-07
## X2_PDRB_Kapita     2.336848e-02 3.916464e-02  0.5966728 5.547975e-01
## X3_Kepadatan      -1.870482e+01 4.953039e+00 -3.7764340 6.316407e-04
## X4_Densitas_Jalan  3.416134e+04 7.881754e+03  4.3342309 1.291062e-04
# =========================================================
# 8. INTERVAL KEPERCAYAAN 95%
# =========================================================

cat("\nINTERVAL KEPERCAYAAN 95%:\n")
## 
## INTERVAL KEPERCAYAAN 95%:
print(confint(model, level = 0.95))
##                           2.5 %        97.5 %
## (Intercept)       -1.762688e+04  1.951876e+03
## X1_Wisnus          1.548252e-04  3.044311e-04
## X2_PDRB_Kapita    -5.631258e-02  1.030495e-01
## X3_Kepadatan      -2.878186e+01 -8.627791e+00
## X4_Densitas_Jalan  1.812579e+04  5.019689e+04
# =========================================================
# 9. UJI t PARSIAL
# =========================================================

hasil_t <- coef(summary(model))

cat("\nUJI t PARSIAL:\n")
## 
## UJI t PARSIAL:
print(hasil_t)
##                        Estimate   Std. Error    t value     Pr(>|t|)
## (Intercept)       -7.837503e+03 4.811652e+03 -1.6288592 1.128536e-01
## X1_Wisnus          2.296282e-04 3.676696e-05  6.2455029 4.692649e-07
## X2_PDRB_Kapita     2.336848e-02 3.916464e-02  0.5966728 5.547975e-01
## X3_Kepadatan      -1.870482e+01 4.953039e+00 -3.7764340 6.316407e-04
## X4_Densitas_Jalan  3.416134e+04 7.881754e+03  4.3342309 1.291062e-04
cat("\nKEPUTUSAN UJI t (ALPHA = 0.05):\n")
## 
## KEPUTUSAN UJI t (ALPHA = 0.05):
p_values <- hasil_t[, 4]

for (i in 2:nrow(hasil_t)) {
  if (p_values[i] < 0.05) {
    cat(rownames(hasil_t)[i], ": SIGNIFIKAN\n")
  } else {
    cat(rownames(hasil_t)[i], ": TIDAK SIGNIFIKAN\n")
  }
}
## X1_Wisnus : SIGNIFIKAN
## X2_PDRB_Kapita : TIDAK SIGNIFIKAN
## X3_Kepadatan : SIGNIFIKAN
## X4_Densitas_Jalan : SIGNIFIKAN
# =========================================================
# 10. UJI F SIMULTAN
# =========================================================

f_stat <- summary(model)$fstatistic

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

p_F <- pf(
  F_value,
  df1,
  df2,
  lower.tail = FALSE
)

cat("\nUJI F SIMULTAN:\n")
## 
## UJI F SIMULTAN:
cat("F hitung =", F_value, "\n")
## F hitung = 20.7243
cat("df1 =", df1, "\n")
## df1 = 4
cat("df2 =", df2, "\n")
## df2 = 33
cat("p-value =", p_F, "\n")
## p-value = 1.275116e-08
if (p_F < 0.05) {
  cat("Keputusan: MODEL SIGNIFIKAN SECARA SIMULTAN\n")
} else {
  cat("Keputusan: MODEL TIDAK SIGNIFIKAN SECARA SIMULTAN\n")
}
## Keputusan: MODEL SIGNIFIKAN SECARA SIMULTAN
# =========================================================
# 11. KOEFISIEN DETERMINASI
# =========================================================

R2 <- summary(model)$r.squared
Adj_R2 <- summary(model)$adj.r.squared
RSE <- summary(model)$sigma

cat("\nKOEFISIEN DETERMINASI:\n")
## 
## KOEFISIEN DETERMINASI:
cat("R-squared =", R2, "\n")
## R-squared = 0.7152649
cat("Adjusted R-squared =", Adj_R2, "\n")
## Adjusted R-squared = 0.6807516
cat("Residual Standard Error =", RSE, "\n")
## Residual Standard Error = 10766.03
cat("\nR-squared (%) =", R2 * 100, "%\n")
## 
## R-squared (%) = 71.52649 %
cat("Adjusted R-squared (%) =", Adj_R2 * 100, "%\n")
## Adjusted R-squared (%) = 68.07516 %
# =========================================================
# 12. RESIDUAL DAN FITTED VALUE
# =========================================================

residuals_model <- residuals(model)
fitted_values <- fitted(model)

cat("\nRESIDUAL:\n")
## 
## RESIDUAL:
print(head(residuals_model))
##          1          2          3          4          5          6 
## -6791.7245 -8316.9470 -7599.5123   379.7865  -208.0417  3514.6146
cat("\nFITTED VALUE:\n")
## 
## FITTED VALUE:
print(head(fitted_values))
##         1         2         3         4         5         6 
## 10344.725 21886.947 13494.512  9234.213  3450.042  5723.385
# =========================================================
# 13. UJI NORMALITAS SHAPIRO-WILK
# =========================================================

shapiro_result <- shapiro.test(residuals_model)

cat("\nUJI SHAPIRO-WILK:\n")
## 
## UJI SHAPIRO-WILK:
print(shapiro_result)
## 
##  Shapiro-Wilk normality test
## 
## data:  residuals_model
## W = 0.78056, p-value = 4.301e-06
# =========================================================
# 14. UJI NORMALITAS TAMBAHAN
# =========================================================

cat("\nUJI ANDERSON-DARLING:\n")
## 
## UJI ANDERSON-DARLING:
print(ad.test(residuals_model))
## 
##  Anderson-Darling normality test
## 
## data:  residuals_model
## A = 1.7705, p-value = 0.0001274
cat("\nUJI LILLIEFORS:\n")
## 
## UJI LILLIEFORS:
print(lillie.test(residuals_model))
## 
##  Lilliefors (Kolmogorov-Smirnov) normality test
## 
## data:  residuals_model
## D = 0.18813, p-value = 0.001592
cat("\nUJI JARQUE-BERA:\n")
## 
## UJI JARQUE-BERA:
print(jarque.test(residuals_model))
## 
##  Jarque-Bera Normality Test
## 
## data:  residuals_model
## JB = 182.83, p-value < 2.2e-16
## alternative hypothesis: greater
# =========================================================
# 15. Q-Q PLOT
# =========================================================

qqnorm(
  residuals_model,
  main = "Normal Q-Q Plot Residual"
)

qqline(residuals_model)

# =========================================================
# 16. UJI HETEROSKEDASTISITAS - BREUSCH-PAGAN
# =========================================================

bp_result <- bptest(model)

cat("\nUJI BREUSCH-PAGAN:\n")
## 
## UJI BREUSCH-PAGAN:
print(bp_result)
## 
##  studentized Breusch-Pagan test
## 
## data:  model
## BP = 14.673, df = 4, p-value = 0.005431
# =========================================================
# 17. UJI GLEJSER
# =========================================================

glejser_model <- lm(
  abs(residuals_model) ~
    X1_Wisnus +
    X2_PDRB_Kapita +
    X3_Kepadatan +
    X4_Densitas_Jalan,
  data = mlr_data
)

cat("\nUJI GLEJSER:\n")
## 
## UJI GLEJSER:
print(summary(glejser_model))
## 
## Call:
## lm(formula = abs(residuals_model) ~ X1_Wisnus + X2_PDRB_Kapita + 
##     X3_Kepadatan + X4_Densitas_Jalan, data = mlr_data)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -10450.6  -2715.5   -737.4   1148.1  23095.1 
## 
## Coefficients:
##                     Estimate Std. Error t value Pr(>|t|)    
## (Intercept)       -6.239e+01  2.613e+03  -0.024   0.9811    
## X1_Wisnus          5.403e-05  1.997e-05   2.706   0.0107 *  
## X2_PDRB_Kapita    -1.066e-02  2.127e-02  -0.501   0.6196    
## X3_Kepadatan      -1.258e+01  2.690e+00  -4.676 4.77e-05 ***
## X4_Densitas_Jalan  2.085e+04  4.280e+03   4.871 2.69e-05 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 5846 on 33 degrees of freedom
## Multiple R-squared:  0.5059, Adjusted R-squared:  0.446 
## F-statistic: 8.448 on 4 and 33 DF,  p-value: 8.285e-05
# =========================================================
# 18. UJI AUTOKORELASI - DURBIN-WATSON
# =========================================================

dw_result <- dwtest(model)

cat("\nUJI DURBIN-WATSON:\n")
## 
## UJI DURBIN-WATSON:
print(dw_result)
## 
##  Durbin-Watson test
## 
## data:  model
## DW = 1.4141, p-value = 0.02288
## alternative hypothesis: true autocorrelation is greater than 0
# =========================================================
# 19. UJI AUTOKORELASI - BREUSCH-GODFREY
# =========================================================

bg_result <- bgtest(
  model,
  order = 1
)

cat("\nUJI BREUSCH-GODFREY:\n")
## 
## UJI BREUSCH-GODFREY:
print(bg_result)
## 
##  Breusch-Godfrey test for serial correlation of order up to 1
## 
## data:  model
## LM test = 3.3211, df = 1, p-value = 0.0684
# =========================================================
# 20. UJI MULTIKOLINEARITAS - VIF
# =========================================================

vif_result <- vif(model)

cat("\nUJI VIF:\n")
## 
## UJI VIF:
print(vif_result)
##         X1_Wisnus    X2_PDRB_Kapita      X3_Kepadatan X4_Densitas_Jalan 
##          1.185370          2.076602         53.277531         48.645249
# =========================================================
# 21. KORELASI X3 DAN X4
# =========================================================

cor_X3_X4 <- cor(
  mlr_data$X3_Kepadatan,
  mlr_data$X4_Densitas_Jalan
)

cat("\nKORELASI X3 DAN X4:\n")
## 
## KORELASI X3 DAN X4:
cat("Korelasi =", cor_X3_X4, "\n")
## Korelasi = 0.9891195
# =========================================================
# 22. UJI RAMSEY RESET
# =========================================================

reset_result <- resettest(
  model,
  power = 2:3,
  type = "fitted"
)

cat("\nUJI RAMSEY RESET:\n")
## 
## UJI RAMSEY RESET:
print(reset_result)
## 
##  RESET test
## 
## data:  model
## RESET = 6.1497, df1 = 2, df2 = 31, p-value = 0.005632
# =========================================================
# 23. COOK'S DISTANCE
# =========================================================

cook_values <- cooks.distance(model)

max_cook <- max(cook_values)
index_cook <- which.max(cook_values)

cat("\nCOOK'S DISTANCE:\n")
## 
## COOK'S DISTANCE:
cat("Nilai maksimum =", max_cook, "\n")
## Nilai maksimum = 115.9722
cat("Indeks observasi =", index_cook, "\n")
## Indeks observasi = 11
plot(
  cook_values,
  type = "h",
  main = "Cook's Distance",
  xlab = "Observasi",
  ylab = "Cook's Distance"
)

abline(
  h = 4 / nrow(mlr_data),
  lty = 2
)

# =========================================================
# 24. LEVERAGE
# =========================================================

leverage <- hatvalues(model)

max_leverage <- max(leverage)
index_leverage <- which.max(leverage)

cat("\nLEVERAGE:\n")
## 
## LEVERAGE:
cat("Nilai maksimum =", max_leverage, "\n")
## Nilai maksimum = 0.9914419
cat("Indeks observasi =", index_leverage, "\n")
## Indeks observasi = 11
# =========================================================
# 25. STUDENTIZED RESIDUAL
# =========================================================

student_res <- rstudent(model)

max_student <- max(student_res)
index_max_student <- which.max(student_res)

min_student <- min(student_res)
index_min_student <- which.min(student_res)

cat("\nSTUDENTIZED RESIDUAL:\n")
## 
## STUDENTIZED RESIDUAL:
cat("Maksimum =", max_student, "\n")
## Maksimum = 14.80753
cat("Indeks maksimum =", index_max_student, "\n")
## Indeks maksimum = 17
cat("Minimum =", min_student, "\n")
## Minimum = -2.39196
cat("Indeks minimum =", index_min_student, "\n")
## Indeks minimum = 11
# =========================================================
# 26. DFBETAS
# =========================================================

dfbetas_model <- dfbetas(model)

cat("\nNILAI MAKSIMUM ABSOLUT DFBETAS:\n")
## 
## NILAI MAKSIMUM ABSOLUT DFBETAS:
print(
  apply(
    abs(dfbetas_model),
    2,
    max
  )
)
##       (Intercept)         X1_Wisnus    X2_PDRB_Kapita      X3_Kepadatan 
##          4.511382          2.690176          1.069265         10.188235 
## X4_Densitas_Jalan 
##         10.633010
# =========================================================
# 27. PLOT DIAGNOSTIK MODEL
# =========================================================

par(mfrow = c(2, 2))

plot(model)
## Warning in sqrt(crit * p * (1 - hh)/hh): NaNs produced
## Warning in sqrt(crit * p * (1 - hh)/hh): NaNs produced

par(mfrow = c(1, 1))


# =========================================================
# 28. STEPWISE REGRESSION
# =========================================================

step_model <- stepAIC(
  model,
  direction = "both",
  trace = TRUE
)
## 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("\nHASIL STEPWISE REGRESSION:\n")
## 
## HASIL STEPWISE REGRESSION:
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
# =========================================================
# 29. PERBANDINGAN MODEL AWAL DAN STEPWISE
# =========================================================

cat("\n================ MODEL AWAL ================\n")
## 
## ================ MODEL AWAL ================
cat("R-squared =", summary(model)$r.squared, "\n")
## R-squared = 0.7152649
cat("Adjusted R-squared =", summary(model)$adj.r.squared, "\n")
## Adjusted R-squared = 0.6807516
cat("AIC =", AIC(model), "\n")
## AIC = 820.0738
cat("\n================ MODEL STEPWISE ================\n")
## 
## ================ MODEL STEPWISE ================
cat("R-squared =", summary(step_model)$r.squared, "\n")
## R-squared = 0.7121931
cat("Adjusted R-squared =", summary(step_model)$adj.r.squared, "\n")
## Adjusted R-squared = 0.6867983
cat("AIC =", AIC(step_model), "\n")
## AIC = 818.4816
# =========================================================
# 30. BOX-COX
# =========================================================

boxcox(
  model,
  main = "Box-Cox Transformation"
)
## Warning: In lm.fit(x, y, offset = offset, singular.ok = singular.ok, ...) :
##  extra argument 'main' will be disregarded

# =========================================================
# 31. RINGKASAN AKHIR
# =========================================================

cat("\n")
cat("=========================================================\n")
## =========================================================
cat("RINGKASAN HASIL ANALISIS MLR\n")
## RINGKASAN HASIL ANALISIS MLR
cat("=========================================================\n")
## =========================================================
cat("Jumlah observasi =", nrow(mlr_data), "\n")
## Jumlah observasi = 38
cat("\nR-squared =", round(R2, 4), "\n")
## 
## R-squared = 0.7153
cat("Adjusted R-squared =", round(Adj_R2, 4), "\n")
## Adjusted R-squared = 0.6808
cat("Residual Standard Error =", round(RSE, 4), "\n")
## Residual Standard Error = 10766.03
cat("\nF-statistic =", round(F_value, 4), "\n")
## 
## F-statistic = 20.7243
cat("p-value F =", p_F, "\n")
## p-value F = 1.275116e-08
cat("\nVIF:\n")
## 
## VIF:
print(round(vif_result, 3))
##         X1_Wisnus    X2_PDRB_Kapita      X3_Kepadatan X4_Densitas_Jalan 
##             1.185             2.077            53.278            48.645
cat("\nKorelasi X3-X4 =", round(cor_X3_X4, 4), "\n")
## 
## Korelasi X3-X4 = 0.9891
cat("\nCook's Distance maksimum =", round(max_cook, 4), "\n")
## 
## Cook's Distance maksimum = 115.9722
cat("Indeks observasi Cook maksimum =", index_cook, "\n")
## Indeks observasi Cook maksimum = 11
cat("\nLeverage maksimum =", round(max_leverage, 4), "\n")
## 
## Leverage maksimum = 0.9914
cat("Indeks observasi leverage maksimum =", index_leverage, "\n")
## Indeks observasi leverage maksimum = 11
cat("\nStudentized residual maksimum =", round(max_student, 4), "\n")
## 
## Studentized residual maksimum = 14.8075
cat("Indeks maksimum =", index_max_student, "\n")
## Indeks maksimum = 17
cat("\nStudentized residual minimum =", round(min_student, 4), "\n")
## 
## Studentized residual minimum = -2.392
cat("Indeks minimum =", index_min_student, "\n")
## Indeks minimum = 11
cat("\n=========================================================\n")
## 
## =========================================================
cat("ANALISIS SELESAI\n")
## ANALISIS SELESAI
cat("=========================================================\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.