# ============================================================
# ANALISIS MULTIPLE LINEAR REGRESSION (MLR)
# DATASET : STUDENT LIFESTYLE DATASET (KAGGLE)
# JUMLAH OBSERVASI : 2000
# VARIABEL DEPENDEN : GPA
# ============================================================
# ============================================================
# 0. PERSIAPAN: LOAD LIBRARY DAN DATA
# ============================================================
library(car) # VIF
## 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.2
library(lmtest) # Breusch-Pagan, Durbin-Watson, Breusch-Godfrey, RESET
## 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.2
##
## Attaching package: 'zoo'
## The following objects are masked from 'package:base':
##
## as.Date, as.Date.numeric
library(nortest) # Anderson-Darling
## Warning: package 'nortest' was built under R version 4.5.2
library(corrplot) # Matriks korelasi
## Warning: package 'corrplot' was built under R version 4.5.3
## corrplot 0.95 loaded
data <- read.csv("C:/Users/user/Downloads/student_lifestyle_dataset.csv")
# Lihat struktur data
str(data)
## 'data.frame': 2000 obs. of 8 variables:
## $ Student_ID : int 1 2 3 4 5 6 7 8 9 10 ...
## $ Study_Hours_Per_Day : num 6.9 5.3 5.1 6.5 8.1 6 8 8.4 5.2 7.7 ...
## $ Extracurricular_Hours_Per_Day : num 3.8 3.5 3.9 2.1 0.6 2.1 0.7 1.8 3.6 0.7 ...
## $ Sleep_Hours_Per_Day : num 8.7 8 9.2 7.2 6.5 8 5.3 5.6 6.3 9.8 ...
## $ Social_Hours_Per_Day : num 2.8 4.2 1.2 1.7 2.2 0.3 5.7 3 4 4.5 ...
## $ Physical_Activity_Hours_Per_Day: num 1.8 3 4.6 6.5 6.6 7.6 4.3 5.2 4.9 1.3 ...
## $ GPA : num 2.99 2.75 2.67 2.88 3.51 2.85 3.08 3.2 2.82 2.76 ...
## $ Stress_Level : chr "Moderate" "Low" "Low" "Moderate" ...
# 6 baris pertama
head(data)
## Student_ID Study_Hours_Per_Day Extracurricular_Hours_Per_Day
## 1 1 6.9 3.8
## 2 2 5.3 3.5
## 3 3 5.1 3.9
## 4 4 6.5 2.1
## 5 5 8.1 0.6
## 6 6 6.0 2.1
## Sleep_Hours_Per_Day Social_Hours_Per_Day Physical_Activity_Hours_Per_Day GPA
## 1 8.7 2.8 1.8 2.99
## 2 8.0 4.2 3.0 2.75
## 3 9.2 1.2 4.6 2.67
## 4 7.2 1.7 6.5 2.88
## 5 6.5 2.2 6.6 3.51
## 6 8.0 0.3 7.6 2.85
## Stress_Level
## 1 Moderate
## 2 Low
## 3 Low
## 4 Moderate
## 5 High
## 6 Moderate
# Jumlah observasi dan variabel
cat("\nJumlah observasi:", nrow(data), "\n")
##
## Jumlah observasi: 2000
cat("Jumlah variabel:", ncol(data), "\n")
## Jumlah variabel: 8
# Keterangan variabel:
# - Student_ID : identitas observasi, tidak digunakan sebagai prediktor
# - Stress_Level: variabel kategorik, tidak digunakan dalam model MLR
# - GPA : variabel dependen (Y)
# ============================================================
# 1. PEMERIKSAAN DATA
# ============================================================
cat("\n")
cat(rep("=", 60), "\n")
## = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = =
cat("1. PEMERIKSAAN DATA\n")
## 1. PEMERIKSAAN DATA
cat(rep("=", 60), "\n")
## = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = =
# Cek missing value
cat("\nJumlah missing value setiap variabel:\n")
##
## Jumlah missing value setiap variabel:
print(colSums(is.na(data)))
## Student_ID Study_Hours_Per_Day
## 0 0
## Extracurricular_Hours_Per_Day Sleep_Hours_Per_Day
## 0 0
## Social_Hours_Per_Day Physical_Activity_Hours_Per_Day
## 0 0
## GPA Stress_Level
## 0 0
# Statistik deskriptif
cat("\nStatistik deskriptif:\n")
##
## Statistik deskriptif:
print(summary(data))
## Student_ID Study_Hours_Per_Day Extracurricular_Hours_Per_Day
## Min. : 1.0 Min. : 5.000 Min. :0.00
## 1st Qu.: 500.8 1st Qu.: 6.300 1st Qu.:1.00
## Median :1000.5 Median : 7.400 Median :2.00
## Mean :1000.5 Mean : 7.476 Mean :1.99
## 3rd Qu.:1500.2 3rd Qu.: 8.700 3rd Qu.:3.00
## Max. :2000.0 Max. :10.000 Max. :4.00
## Sleep_Hours_Per_Day Social_Hours_Per_Day Physical_Activity_Hours_Per_Day
## Min. : 5.000 Min. :0.000 Min. : 0.000
## 1st Qu.: 6.200 1st Qu.:1.200 1st Qu.: 2.400
## Median : 7.500 Median :2.600 Median : 4.100
## Mean : 7.501 Mean :2.705 Mean : 4.328
## 3rd Qu.: 8.800 3rd Qu.:4.100 3rd Qu.: 6.100
## Max. :10.000 Max. :6.000 Max. :13.000
## GPA Stress_Level
## Min. :2.240 Length:2000
## 1st Qu.:2.900 Class :character
## Median :3.110 Mode :character
## Mean :3.116
## 3rd Qu.:3.330
## Max. :4.000
# ============================================================
# 2. PEMERIKSAAN TOTAL WAKTU HARIAN
# ============================================================
cat("\n")
cat(rep("=", 60), "\n")
## = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = =
cat("2. PEMERIKSAAN TOTAL WAKTU HARIAN\n")
## 2. PEMERIKSAAN TOTAL WAKTU HARIAN
cat(rep("=", 60), "\n")
## = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = =
# Lima aktivitas harian
data$Total_Hours <- rowSums(data[, c("Study_Hours_Per_Day", "Extracurricular_Hours_Per_Day", "Sleep_Hours_Per_Day", "Social_Hours_Per_Day", "Physical_Activity_Hours_Per_Day")])
cat("\nRingkasan total jam harian:\n")
##
## Ringkasan total jam harian:
print(summary(data$Total_Hours))
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 24 24 24 24 24 24
cat("\nApakah seluruh observasi berjumlah 24 jam?\n")
##
## Apakah seluruh observasi berjumlah 24 jam?
print(all(abs(data$Total_Hours - 24) < 1e-10))
## [1] TRUE
# Karena kelima aktivitas berjumlah tepat 24 jam:
#
# Study + Extracurricular + Sleep + Social + Physical Activity = 24
#
# maka kelima variabel tidak dapat dimasukkan sekaligus ke dalam
# model dengan intercept (perfect multicollinearity).
#
# Physical Activity dijadikan variabel referensi, sehingga model
# menggunakan empat variabel:
# - Study Hours
# - Extracurricular Hours
# - Sleep Hours
# - Social Hours
# Pembuktian perfect multicollinearity
model_check <- lm(GPA ~ Study_Hours_Per_Day + Extracurricular_Hours_Per_Day + Sleep_Hours_Per_Day + Social_Hours_Per_Day + Physical_Activity_Hours_Per_Day, data = data)
cat("\nKoefisien jika kelima variabel dimasukkan sekaligus:\n")
##
## Koefisien jika kelima variabel dimasukkan sekaligus:
print(coef(model_check))
## (Intercept) Study_Hours_Per_Day
## 2.007303084 0.154384498
## Extracurricular_Hours_Per_Day Sleep_Hours_Per_Day
## -0.007496213 -0.004548773
## Social_Hours_Per_Day Physical_Activity_Hours_Per_Day
## 0.001312227 NA
cat("(Koefisien NA menunjukkan adanya perfect multicollinearity)\n")
## (Koefisien NA menunjukkan adanya perfect multicollinearity)
# ============================================================
# 3. STATISTIK DESKRIPTIF DAN EKSPLORASI
# ============================================================
cat("\n")
cat(rep("=", 60), "\n")
## = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = =
cat("3. STATISTIK DESKRIPTIF DAN EKSPLORASI\n")
## 3. STATISTIK DESKRIPTIF DAN EKSPLORASI
cat(rep("=", 60), "\n")
## = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = =
# Variabel yang digunakan dalam analisis
analysis_data <- data[, c("Study_Hours_Per_Day", "Extracurricular_Hours_Per_Day", "Sleep_Hours_Per_Day", "Social_Hours_Per_Day", "GPA")]
# Statistik deskriptif
print(summary(analysis_data))
## Study_Hours_Per_Day Extracurricular_Hours_Per_Day Sleep_Hours_Per_Day
## Min. : 5.000 Min. :0.00 Min. : 5.000
## 1st Qu.: 6.300 1st Qu.:1.00 1st Qu.: 6.200
## Median : 7.400 Median :2.00 Median : 7.500
## Mean : 7.476 Mean :1.99 Mean : 7.501
## 3rd Qu.: 8.700 3rd Qu.:3.00 3rd Qu.: 8.800
## Max. :10.000 Max. :4.00 Max. :10.000
## Social_Hours_Per_Day GPA
## Min. :0.000 Min. :2.240
## 1st Qu.:1.200 1st Qu.:2.900
## Median :2.600 Median :3.110
## Mean :2.705 Mean :3.116
## 3rd Qu.:4.100 3rd Qu.:3.330
## Max. :6.000 Max. :4.000
# Matriks korelasi
cat("\nMatriks Korelasi:\n")
##
## Matriks Korelasi:
cor_matrix <- cor(analysis_data)
print(round(cor_matrix, 3))
## Study_Hours_Per_Day Extracurricular_Hours_Per_Day
## Study_Hours_Per_Day 1.000 -0.003
## Extracurricular_Hours_Per_Day -0.003 1.000
## Sleep_Hours_Per_Day 0.027 0.009
## Social_Hours_Per_Day -0.138 -0.139
## GPA 0.734 -0.032
## Sleep_Hours_Per_Day Social_Hours_Per_Day GPA
## Study_Hours_Per_Day 0.027 -0.138 0.734
## Extracurricular_Hours_Per_Day 0.009 -0.139 -0.032
## Sleep_Hours_Per_Day 1.000 -0.194 -0.004
## Social_Hours_Per_Day -0.194 1.000 -0.086
## GPA -0.004 -0.086 1.000
# Visualisasi 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(analysis_data, main = "Scatter Plot Matrix", pch = 19)
# ============================================================
# 4. ESTIMASI MODEL AWAL (FULL MODEL)
# ============================================================
cat("\n")
cat(rep("=", 60), "\n")
## = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = =
cat("4. ESTIMASI MODEL AWAL (FULL MODEL)\n")
## 4. ESTIMASI MODEL AWAL (FULL MODEL)
cat(rep("=", 60), "\n")
## = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = =
# Physical Activity tidak dimasukkan (variabel referensi).
model <- lm(GPA ~ Study_Hours_Per_Day + Extracurricular_Hours_Per_Day + Sleep_Hours_Per_Day + Social_Hours_Per_Day, data = data)
# Ringkasan model
cat("\nRingkasan Model Awal:\n")
##
## Ringkasan Model Awal:
print(summary(model))
##
## Call:
## lm(formula = GPA ~ Study_Hours_Per_Day + Extracurricular_Hours_Per_Day +
## Sleep_Hours_Per_Day + Social_Hours_Per_Day, data = data)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.59447 -0.13516 -0.00293 0.13439 0.78728
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 2.007303 0.037815 53.082 <2e-16 ***
## Study_Hours_Per_Day 0.154384 0.003213 48.046 <2e-16 ***
## Extracurricular_Hours_Per_Day -0.007496 0.003960 -1.893 0.0585 .
## Sleep_Hours_Per_Day -0.004549 0.003161 -1.439 0.1503
## Social_Hours_Per_Day 0.001312 0.002788 0.471 0.6380
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.2026 on 1995 degrees of freedom
## Multiple R-squared: 0.541, Adjusted R-squared: 0.54
## F-statistic: 587.8 on 4 and 1995 DF, p-value: < 2.2e-16
# Koefisien regresi
coef_model <- summary(model)$coefficients
coef_df <- data.frame(
Variabel = rownames(coef_model),
Koefisien = round(coef_model[, 1], 4),
Std_Error = round(coef_model[, 2], 4),
t_value = round(coef_model[, 3], 4),
p_value = round(coef_model[, 4], 4)
)
cat("\nKoefisien Regresi Model Awal:\n")
##
## Koefisien Regresi Model Awal:
print(coef_df)
## Variabel Koefisien Std_Error
## (Intercept) (Intercept) 2.0073 0.0378
## Study_Hours_Per_Day Study_Hours_Per_Day 0.1544 0.0032
## Extracurricular_Hours_Per_Day Extracurricular_Hours_Per_Day -0.0075 0.0040
## Sleep_Hours_Per_Day Sleep_Hours_Per_Day -0.0045 0.0032
## Social_Hours_Per_Day Social_Hours_Per_Day 0.0013 0.0028
## t_value p_value
## (Intercept) 53.0818 0.0000
## Study_Hours_Per_Day 48.0464 0.0000
## Extracurricular_Hours_Per_Day -1.8931 0.0585
## Sleep_Hours_Per_Day -1.4389 0.1503
## Social_Hours_Per_Day 0.4706 0.6380
# Interval kepercayaan 95%
cat("\nInterval Kepercayaan 95%:\n")
##
## Interval Kepercayaan 95%:
print(round(confint(model, level = 0.95), 4))
## 2.5 % 97.5 %
## (Intercept) 1.9331 2.0815
## Study_Hours_Per_Day 0.1481 0.1607
## Extracurricular_Hours_Per_Day -0.0153 0.0003
## Sleep_Hours_Per_Day -0.0107 0.0017
## Social_Hours_Per_Day -0.0042 0.0068
# ============================================================
# 5. PEMILIHAN MODEL (FORWARD, BACKWARD, STEPWISE)
# ============================================================
cat("\n")
cat(rep("=", 60), "\n")
## = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = =
cat("5. PEMILIHAN MODEL\n")
## 5. PEMILIHAN MODEL
cat(rep("=", 60), "\n")
## = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = =
# ------------------------------------------------------------
# Forward
# ------------------------------------------------------------
model_null <- lm(GPA ~ 1, data = data)
model_forward <- step(model_null, scope = formula(model), direction = "forward", trace = 1)
## Start: AIC=-4832.62
## GPA ~ 1
##
## Df Sum of Sq RSS AIC
## + Study_Hours_Per_Day 1 96.195 82.128 -6381.3
## + Social_Hours_Per_Day 1 1.309 177.014 -4845.4
## + Extracurricular_Hours_Per_Day 1 0.185 178.138 -4832.7
## <none> 178.323 -4832.6
## + Sleep_Hours_Per_Day 1 0.003 178.319 -4830.7
##
## Step: AIC=-6381.26
## GPA ~ Study_Hours_Per_Day
##
## Df Sum of Sq RSS AIC
## + Extracurricular_Hours_Per_Day 1 0.16310 81.965 -6383.2
## + Sleep_Hours_Per_Day 1 0.10194 82.026 -6381.7
## <none> 82.128 -6381.3
## + Social_Hours_Per_Day 1 0.04394 82.084 -6380.3
##
## Step: AIC=-6383.23
## GPA ~ Study_Hours_Per_Day + Extracurricular_Hours_Per_Day
##
## Df Sum of Sq RSS AIC
## + Sleep_Hours_Per_Day 1 0.099666 81.865 -6383.7
## <none> 81.965 -6383.2
## + Social_Hours_Per_Day 1 0.023808 81.941 -6381.8
##
## Step: AIC=-6383.66
## GPA ~ Study_Hours_Per_Day + Extracurricular_Hours_Per_Day + Sleep_Hours_Per_Day
##
## Df Sum of Sq RSS AIC
## <none> 81.865 -6383.7
## + Social_Hours_Per_Day 1 0.0090872 81.856 -6381.9
cat("\n--- MODEL FORWARD ---\n")
##
## --- MODEL FORWARD ---
print(formula(model_forward))
## GPA ~ Study_Hours_Per_Day + Extracurricular_Hours_Per_Day + Sleep_Hours_Per_Day
# ------------------------------------------------------------
# Backward
# ------------------------------------------------------------
model_backward <- step(model, direction = "backward", trace = 1)
## Start: AIC=-6381.89
## GPA ~ Study_Hours_Per_Day + Extracurricular_Hours_Per_Day + Sleep_Hours_Per_Day +
## Social_Hours_Per_Day
##
## Df Sum of Sq RSS AIC
## - Social_Hours_Per_Day 1 0.009 81.865 -6383.7
## <none> 81.856 -6381.9
## - Sleep_Hours_Per_Day 1 0.085 81.941 -6381.8
## - Extracurricular_Hours_Per_Day 1 0.147 82.003 -6380.3
## - Study_Hours_Per_Day 1 94.717 176.573 -4846.3
##
## Step: AIC=-6383.66
## GPA ~ Study_Hours_Per_Day + Extracurricular_Hours_Per_Day + Sleep_Hours_Per_Day
##
## Df Sum of Sq RSS AIC
## <none> 81.865 -6383.7
## - Sleep_Hours_Per_Day 1 0.100 81.965 -6383.2
## - Extracurricular_Hours_Per_Day 1 0.161 82.026 -6381.7
## - Study_Hours_Per_Day 1 96.270 178.135 -4830.7
cat("\n--- MODEL BACKWARD ---\n")
##
## --- MODEL BACKWARD ---
print(formula(model_backward))
## GPA ~ Study_Hours_Per_Day + Extracurricular_Hours_Per_Day + Sleep_Hours_Per_Day
# ------------------------------------------------------------
# Stepwise / Both
# ------------------------------------------------------------
model_stepwise <- step(model, direction = "both", trace = 1)
## Start: AIC=-6381.89
## GPA ~ Study_Hours_Per_Day + Extracurricular_Hours_Per_Day + Sleep_Hours_Per_Day +
## Social_Hours_Per_Day
##
## Df Sum of Sq RSS AIC
## - Social_Hours_Per_Day 1 0.009 81.865 -6383.7
## <none> 81.856 -6381.9
## - Sleep_Hours_Per_Day 1 0.085 81.941 -6381.8
## - Extracurricular_Hours_Per_Day 1 0.147 82.003 -6380.3
## - Study_Hours_Per_Day 1 94.717 176.573 -4846.3
##
## Step: AIC=-6383.66
## GPA ~ Study_Hours_Per_Day + Extracurricular_Hours_Per_Day + Sleep_Hours_Per_Day
##
## Df Sum of Sq RSS AIC
## <none> 81.865 -6383.7
## - Sleep_Hours_Per_Day 1 0.100 81.965 -6383.2
## + Social_Hours_Per_Day 1 0.009 81.856 -6381.9
## - Extracurricular_Hours_Per_Day 1 0.161 82.026 -6381.7
## - Study_Hours_Per_Day 1 96.270 178.135 -4830.7
cat("\n--- MODEL STEPWISE / BOTH ---\n")
##
## --- MODEL STEPWISE / BOTH ---
print(formula(model_stepwise))
## GPA ~ Study_Hours_Per_Day + Extracurricular_Hours_Per_Day + Sleep_Hours_Per_Day
# ============================================================
# 6. PERBANDINGAN MODEL DAN MODEL AKHIR
# ============================================================
cat("\n")
cat(rep("=", 60), "\n")
## = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = =
cat("6. PERBANDINGAN MODEL DAN MODEL AKHIR\n")
## 6. PERBANDINGAN MODEL DAN MODEL AKHIR
cat(rep("=", 60), "\n")
## = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = =
# Fungsi AICc (k = jumlah koefisien + 1 parameter varians)
AICc <- function(model) {
n <- nobs(model)
k <- length(coef(model)) + 1
AIC(model) + (2 * k * (k + 1)) / (n - k - 1)
}
model_compare <- data.frame(
Model = c("Model Awal", "Forward", "Backward", "Stepwise"),
Jumlah_X = c(length(coef(model)) - 1, length(coef(model_forward)) - 1, length(coef(model_backward)) - 1, length(coef(model_stepwise)) - 1),
AIC = c(AIC(model), AIC(model_forward), AIC(model_backward), AIC(model_stepwise)),
BIC = c(BIC(model), BIC(model_forward), BIC(model_backward), BIC(model_stepwise)),
AICc = c(AICc(model), AICc(model_forward), AICc(model_backward), AICc(model_stepwise)),
Adjusted_R2 = c(summary(model)$adj.r.squared, summary(model_forward)$adj.r.squared, summary(model_backward)$adj.r.squared, summary(model_stepwise)$adj.r.squared)
)
model_compare[, 3:6] <- round(model_compare[, 3:6], 4)
print(model_compare)
## Model Jumlah_X AIC BIC AICc Adjusted_R2
## 1 Model Awal 4 -704.1322 -670.5267 -704.09 0.5400
## 2 Forward 3 -705.9101 -677.9056 -705.88 0.5402
## 3 Backward 3 -705.9101 -677.9056 -705.88 0.5402
## 4 Stepwise 3 -705.9101 -677.9056 -705.88 0.5402
# ------------------------------------------------------------
# Pemilihan kriteria berdasarkan ukuran data
# ------------------------------------------------------------
# AICc digunakan jika n/k kecil (aturan umum: n/k < 40), karena pada
# sampel kecil AIC cenderung memilih model yang terlalu kompleks.
# Jika n/k besar, AICc hampir sama dengan AIC sehingga AIC digunakan.
n_data <- nrow(data)
k_full <- length(coef(model)) + 1
rasio_n_k <- n_data / k_full
kriteria <- ifelse(rasio_n_k < 40, "AICc", "AIC")
cat("\nn =", n_data, "| k =", k_full, "| n/k =", round(rasio_n_k, 2), "\n")
##
## n = 2000 | k = 6 | n/k = 333.33
cat("Kriteria yang digunakan untuk memilih model akhir:", kriteria, "\n")
## Kriteria yang digunakan untuk memilih model akhir: AIC
if (kriteria == "AICc") {
idx_final <- which.min(model_compare$AICc)
} else {
idx_final <- which.min(model_compare$AIC)
}
# ------------------------------------------------------------
# Model akhir
# ------------------------------------------------------------
daftar_model <- list(model, model_forward, model_backward, model_stepwise)
model_final <- daftar_model[[idx_final]]
cat("\nModel akhir terpilih:", model_compare$Model[idx_final], "\n")
##
## Model akhir terpilih: Forward
cat("\nFormula Model Akhir:\n")
##
## Formula Model Akhir:
print(formula(model_final))
## GPA ~ Study_Hours_Per_Day + Extracurricular_Hours_Per_Day + Sleep_Hours_Per_Day
cat("\nRingkasan Model Akhir:\n")
##
## Ringkasan Model Akhir:
print(summary(model_final))
##
## Call:
## lm(formula = GPA ~ Study_Hours_Per_Day + Extracurricular_Hours_Per_Day +
## Sleep_Hours_Per_Day, data = data)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.59676 -0.13630 -0.00346 0.13466 0.78628
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 2.015074 0.034013 59.244 <2e-16 ***
## Study_Hours_Per_Day 0.154177 0.003182 48.448 <2e-16 ***
## Extracurricular_Hours_Per_Day -0.007760 0.003919 -1.980 0.0478 *
## Sleep_Hours_Per_Day -0.004835 0.003102 -1.559 0.1192
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.2025 on 1996 degrees of freedom
## Multiple R-squared: 0.5409, Adjusted R-squared: 0.5402
## F-statistic: 783.9 on 3 and 1996 DF, p-value: < 2.2e-16
cat("\nInterval Kepercayaan 95% Model Akhir:\n")
##
## Interval Kepercayaan 95% Model Akhir:
print(round(confint(model_final, level = 0.95), 4))
## 2.5 % 97.5 %
## (Intercept) 1.9484 2.0818
## Study_Hours_Per_Day 0.1479 0.1604
## Extracurricular_Hours_Per_Day -0.0154 -0.0001
## Sleep_Hours_Per_Day -0.0109 0.0012
# ============================================================
# 7. UJI HIPOTESIS MODEL AKHIR
# ============================================================
cat("\n")
cat(rep("=", 60), "\n")
## = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = =
cat("7. UJI HIPOTESIS MODEL AKHIR\n")
## 7. UJI HIPOTESIS MODEL AKHIR
cat(rep("=", 60), "\n")
## = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = =
final_summary <- summary(model_final)
# ------------------------------------------------------------
# 7.1 Uji F (simultan)
# ------------------------------------------------------------
f_stat <- final_summary$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("\n--- UJI F (SIMULTAN) ---\n")
##
## --- UJI F (SIMULTAN) ---
cat("H0: seluruh koefisien slope = 0\n")
## H0: seluruh koefisien slope = 0
cat("H1: minimal terdapat satu koefisien slope ≠ 0\n")
## H1: minimal terdapat satu koefisien slope ≠ 0
cat("\nF-statistic:", round(f_value, 4), "\n")
##
## F-statistic: 783.9318
cat("df1:", df1, "\n")
## df1: 3
cat("df2:", df2, "\n")
## df2: 1996
cat("p-value:", format(f_pvalue, scientific = TRUE), "\n")
## p-value: 0e+00
if (f_pvalue < 0.05) {
cat("Kesimpulan: Tolak H0 pada α = 5%. Model signifikan secara simultan.\n")
} else {
cat("Kesimpulan: Gagal menolak H0 pada α = 5%.\n")
}
## Kesimpulan: Tolak H0 pada α = 5%. Model signifikan secara simultan.
# ------------------------------------------------------------
# 7.2 Uji t (parsial)
# ------------------------------------------------------------
tabel_t <- final_summary$coefficients
cat("\n--- UJI t (PARSIAL) ---\n")
##
## --- UJI t (PARSIAL) ---
cat("H0: βj = 0\n")
## H0: βj = 0
cat("H1: βj ≠ 0\n")
## H1: βj ≠ 0
cat("\nTabel koefisien model akhir:\n")
##
## Tabel koefisien model akhir:
print(round(tabel_t, 4))
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 2.0151 0.0340 59.2445 0.0000
## Study_Hours_Per_Day 0.1542 0.0032 48.4482 0.0000
## Extracurricular_Hours_Per_Day -0.0078 0.0039 -1.9802 0.0478
## Sleep_Hours_Per_Day -0.0048 0.0031 -1.5588 0.1192
cat("\nKeputusan (α = 5%):\n")
##
## Keputusan (α = 5%):
for (i in 2:nrow(tabel_t)) {
var_name <- rownames(tabel_t)[i]
t_val <- tabel_t[i, 3]
p_val <- tabel_t[i, 4]
keputusan <- ifelse(p_val < 0.05, "signifikan", "tidak signifikan")
cat(sprintf("%-35s : t = %8.4f, p-value = %8.4f -> %s\n", var_name, t_val, p_val, keputusan))
}
## Study_Hours_Per_Day : t = 48.4482, p-value = 0.0000 -> signifikan
## Extracurricular_Hours_Per_Day : t = -1.9802, p-value = 0.0478 -> signifikan
## Sleep_Hours_Per_Day : t = -1.5588, p-value = 0.1192 -> tidak signifikan
# ------------------------------------------------------------
# 7.3 Koefisien determinasi
# ------------------------------------------------------------
r_squared <- final_summary$r.squared
adj_r_squared <- final_summary$adj.r.squared
resid_se <- final_summary$sigma
cat("\n--- KOEFISIEN DETERMINASI ---\n")
##
## --- KOEFISIEN DETERMINASI ---
cat("R-squared:", round(r_squared, 4), "\n")
## R-squared: 0.5409
cat("Adjusted R-squared:", round(adj_r_squared, 4), "\n")
## Adjusted R-squared: 0.5402
cat("Residual Standard Error:", round(resid_se, 4), "\n")
## Residual Standard Error: 0.2025
cat("Model menjelaskan sekitar", round(r_squared * 100, 2), "% variasi GPA.\n")
## Model menjelaskan sekitar 54.09 % variasi GPA.
# ------------------------------------------------------------
# 7.4 Persamaan regresi
# ------------------------------------------------------------
final_coef <- coef(model_final)
cat("\n--- PERSAMAAN REGRESI MODEL AKHIR ---\n")
##
## --- PERSAMAAN REGRESI MODEL AKHIR ---
cat("GPA =", round(final_coef[1], 4))
## GPA = 2.0151
if (length(final_coef) > 1) {
for (i in 2:length(final_coef)) {
if (final_coef[i] >= 0) {
cat(" + ", round(final_coef[i], 4), "*", names(final_coef)[i], sep = "")
} else {
cat(" - ", round(abs(final_coef[i]), 4), "*", names(final_coef)[i], sep = "")
}
}
}
## + 0.1542*Study_Hours_Per_Day - 0.0078*Extracurricular_Hours_Per_Day - 0.0048*Sleep_Hours_Per_Day
cat("\n")
# Catatan interpretasi:
# Karena Physical Activity menjadi variabel referensi dalam alokasi
# waktu 24 jam, koefisien setiap variabel waktu dibaca sebagai
# perubahan GPA jika 1 jam dialihkan dari Physical Activity ke
# variabel tersebut, dengan variabel lain dalam model dianggap tetap.
# ============================================================
# 8. UJI ASUMSI MODEL AKHIR
# ============================================================
cat("\n")
cat(rep("=", 60), "\n")
## = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = =
cat("8. UJI ASUMSI MODEL AKHIR\n")
## 8. UJI ASUMSI MODEL AKHIR
cat(rep("=", 60), "\n")
## = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = =
# Residual
residuals_model <- residuals(model_final)
# Fitted value
fitted_values <- fitted(model_final)
# Standardized residual
standardized_resid <- rstandard(model_final)
# ------------------------------------------------------------
# 8.1 UJI NORMALITAS RESIDUAL
# ------------------------------------------------------------
cat("\n--- 8.1 UJI NORMALITAS RESIDUAL ---\n")
##
## --- 8.1 UJI NORMALITAS RESIDUAL ---
# Anderson-Darling dipilih karena:
# - jumlah data besar (n = 2000)
# - sensitif terhadap penyimpangan di bagian ekor distribusi
# - tersedia pada package nortest
#
# Shapiro-Wilk lebih sesuai untuk sampel kecil, dan Kolmogorov-Smirnov
# kurang sensitif untuk parameter yang diestimasi dari data.
ad_test <- ad.test(residuals_model)
cat("\nAnderson-Darling Test:\n")
##
## Anderson-Darling Test:
cat("H0: residual berdistribusi normal\n")
## H0: residual berdistribusi normal
cat("H1: residual tidak berdistribusi normal\n")
## H1: residual tidak berdistribusi normal
cat("A =", round(ad_test$statistic, 4), "\n")
## A = 0.4531
cat("p-value =", format(ad_test$p.value, scientific = TRUE), "\n")
## p-value = 2.711342e-01
if (ad_test$p.value > 0.05) {
cat("Kesimpulan: Tidak terdapat bukti yang cukup untuk menyatakan residual tidak normal.\n")
} else {
cat("Kesimpulan: Terdapat bukti bahwa residual tidak berdistribusi normal.\n")
}
## Kesimpulan: Tidak terdapat bukti yang cukup untuk menyatakan residual tidak normal.
# ------------------------------------------------------------
# 8.2 UJI HETEROSKEDASTISITAS
# ------------------------------------------------------------
cat("\n--- 8.2 UJI HETEROSKEDASTISITAS ---\n")
##
## --- 8.2 UJI HETEROSKEDASTISITAS ---
bp_test <- bptest(model_final)
cat("\nBreusch-Pagan Test:\n")
##
## Breusch-Pagan Test:
cat("H0: varians residual homogen\n")
## H0: varians residual homogen
cat("H1: varians residual tidak homogen\n")
## H1: varians residual tidak homogen
cat("BP =", round(bp_test$statistic, 4), "\n")
## BP = 1.9558
cat("p-value =", round(bp_test$p.value, 4), "\n")
## p-value = 0.5816
if (bp_test$p.value > 0.05) {
cat("Kesimpulan: Tidak terdapat bukti heteroskedastisitas.\n")
} else {
cat("Kesimpulan: Terdapat bukti heteroskedastisitas.\n")
}
## Kesimpulan: Tidak terdapat bukti heteroskedastisitas.
# ------------------------------------------------------------
# 8.3 UJI AUTOKORELASI (PELENGKAP)
# ------------------------------------------------------------
cat("\n--- 8.3 UJI AUTOKORELASI (PELENGKAP) ---\n")
##
## --- 8.3 UJI AUTOKORELASI (PELENGKAP) ---
# Data merupakan data cross-sectional mahasiswa, sehingga tidak ada
# urutan waktu yang bermakna dan Student_ID hanya identitas observasi.
# Independensi observasi terutama dinilai dari desain pengumpulan data.
# Durbin-Watson dan Breusch-Godfrey ditampilkan sebagai pelengkap saja.
dw_test <- dwtest(model_final)
cat("\nDurbin-Watson Test:\n")
##
## Durbin-Watson Test:
cat("DW =", round(dw_test$statistic, 4), "\n")
## DW = 2.0058
cat("p-value =", round(dw_test$p.value, 4), "\n")
## p-value = 0.5514
bg_test <- bgtest(model_final, order = 1)
cat("\nBreusch-Godfrey Test (Lag 1):\n")
##
## Breusch-Godfrey Test (Lag 1):
cat("LM =", round(bg_test$statistic, 4), "\n")
## LM = 0.0205
cat("p-value =", round(bg_test$p.value, 4), "\n")
## p-value = 0.8862
if (bg_test$p.value > 0.05) {
cat("Kesimpulan: Tidak terdapat bukti autokorelasi berdasarkan urutan data.\n")
} else {
cat("Kesimpulan: Terdapat indikasi autokorelasi berdasarkan urutan data.\n")
}
## Kesimpulan: Tidak terdapat bukti autokorelasi berdasarkan urutan data.
# ------------------------------------------------------------
# 8.4 UJI MULTIKOLINEARITAS
# ------------------------------------------------------------
cat("\n--- 8.4 UJI MULTIKOLINEARITAS ---\n")
##
## --- 8.4 UJI MULTIKOLINEARITAS ---
if (length(coef(model_final)) > 2) {
vif_values <- vif(model_final)
cat("\nVIF Values:\n")
print(round(vif_values, 4))
cat("\nInterpretasi VIF:\n")
for (i in seq_along(vif_values)) {
vif_val <- vif_values[i]
var_name <- names(vif_values)[i]
status <- ifelse(vif_val < 5, "Tidak menunjukkan multikolinearitas", ifelse(vif_val < 10, "Multikolinearitas moderat", "Multikolinearitas serius"))
cat(sprintf("%-35s : VIF = %7.4f -> %s\n", var_name, vif_val, status))
}
# Tolerance
tolerance <- 1 / vif_values
cat("\nTolerance Values:\n")
print(round(tolerance, 4))
vif_maks <- round(max(vif_values), 4)
} else {
cat("\nModel akhir hanya memiliki satu prediktor, sehingga VIF tidak dihitung.\n")
vif_maks <- NA
}
##
## VIF Values:
## Study_Hours_Per_Day Extracurricular_Hours_Per_Day
## 1.0007 1.0001
## Sleep_Hours_Per_Day
## 1.0008
##
## Interpretasi VIF:
## Study_Hours_Per_Day : VIF = 1.0007 -> Tidak menunjukkan multikolinearitas
## Extracurricular_Hours_Per_Day : VIF = 1.0001 -> Tidak menunjukkan multikolinearitas
## Sleep_Hours_Per_Day : VIF = 1.0008 -> Tidak menunjukkan multikolinearitas
##
## Tolerance Values:
## Study_Hours_Per_Day Extracurricular_Hours_Per_Day
## 0.9993 0.9999
## Sleep_Hours_Per_Day
## 0.9992
# ------------------------------------------------------------
# 8.5 UJI LINEARITAS / SPESIFIKASI MODEL
# ------------------------------------------------------------
cat("\n--- 8.5 UJI LINEARITAS / SPESIFIKASI MODEL ---\n")
##
## --- 8.5 UJI LINEARITAS / SPESIFIKASI MODEL ---
reset_test <- resettest(model_final, power = 2:3, type = "fitted")
cat("\nRamsey RESET Test:\n")
##
## Ramsey RESET Test:
cat("H0: spesifikasi model linear sudah tepat\n")
## H0: spesifikasi model linear sudah tepat
cat("H1: terdapat kesalahan spesifikasi model\n")
## H1: terdapat kesalahan spesifikasi model
cat("F =", round(reset_test$statistic, 4), "\n")
## F = 0.4191
cat("p-value =", round(reset_test$p.value, 4), "\n")
## p-value = 0.6577
if (reset_test$p.value > 0.05) {
cat("Kesimpulan: Tidak terdapat bukti adanya kesalahan spesifikasi non-linear.\n")
} else {
cat("Kesimpulan: Terdapat indikasi masalah spesifikasi model.\n")
}
## Kesimpulan: Tidak terdapat bukti adanya kesalahan spesifikasi non-linear.
# Catatan:
# Jika ada uji asumsi yang tidak terpenuhi, data tidak diperbaiki
# (tanpa transformasi).
# ============================================================
# 9. DIAGNOSTIK OUTLIER DAN INFLUENTIAL POINTS
# ============================================================
cat("\n")
cat(rep("=", 60), "\n")
## = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = =
cat("9. DIAGNOSTIK OUTLIER DAN INFLUENTIAL POINTS\n")
## 9. DIAGNOSTIK OUTLIER DAN INFLUENTIAL POINTS
cat(rep("=", 60), "\n")
## = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = =
n_obs <- nrow(data)
k_coef <- length(coef(model_final))
# ------------------------------------------------------------
# 9.1 Cook's Distance
# ------------------------------------------------------------
cooks_d <- cooks.distance(model_final)
cat("\nCook's Distance:\n")
##
## Cook's Distance:
cat("Nilai maksimum:", round(max(cooks_d), 4), "\n")
## Nilai maksimum: 0.0115
cat("Jumlah observasi dengan Cook's D > 1:", sum(cooks_d > 1), "\n")
## Jumlah observasi dengan Cook's D > 1: 0
# ------------------------------------------------------------
# 9.2 Leverage
# ------------------------------------------------------------
leverage <- hatvalues(model_final)
leverage_threshold <- 2 * k_coef / n_obs
cat("\nLeverage (Hat Values):\n")
##
## Leverage (Hat Values):
cat("Nilai maksimum:", round(max(leverage), 4), "\n")
## Nilai maksimum: 0.0047
cat("Threshold 2(k+1)/n:", round(leverage_threshold, 4), "\n")
## Threshold 2(k+1)/n: 0.004
cat("Jumlah observasi dengan leverage > threshold:", sum(leverage > leverage_threshold), "\n")
## Jumlah observasi dengan leverage > threshold: 19
# ------------------------------------------------------------
# 9.3 Studentized Residual
# ------------------------------------------------------------
std_resid <- rstudent(model_final)
cat("\nStudentized Residuals:\n")
##
## Studentized Residuals:
cat("Nilai maksimum absolut:", round(max(abs(std_resid)), 4), "\n")
## Nilai maksimum absolut: 3.9022
cat("Jumlah observasi dengan |studentized residual| > 3:", sum(abs(std_resid) > 3), "\n")
## Jumlah observasi dengan |studentized residual| > 3: 3
# ------------------------------------------------------------
# 9.4 DFBETAS
# ------------------------------------------------------------
dfbetas_val <- dfbetas(model_final)
# - |DFBETAS| > 2/sqrt(n): batas yang disesuaikan dengan ukuran sampel
dfbetas_threshold <- 2 / sqrt(n_obs)
cat("\nDFBETAS:\n")
##
## DFBETAS:
cat("Jumlah observasi dengan minimal satu |DFBETAS| > 1:", sum(apply(abs(dfbetas_val) > 1, 1, any)), "\n")
## Jumlah observasi dengan minimal satu |DFBETAS| > 1: 0
cat("Threshold 2/sqrt(n):", round(dfbetas_threshold, 4), "\n")
## Threshold 2/sqrt(n): 0.0447
cat("Jumlah observasi dengan minimal satu |DFBETAS| > 2/sqrt(n):", sum(apply(abs(dfbetas_val) > dfbetas_threshold, 1, any)), "\n")
## Jumlah observasi dengan minimal satu |DFBETAS| > 2/sqrt(n): 304
cat("\nCatatan: observasi yang terdeteksi tidak otomatis dihapus.\n")
##
## Catatan: observasi yang terdeteksi tidak otomatis dihapus.
# ============================================================
# 10. VISUALISASI DIAGNOSTIK
# ============================================================
cat("\n")
cat(rep("=", 60), "\n")
## = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = =
cat("10. VISUALISASI DIAGNOSTIK\n")
## 10. VISUALISASI DIAGNOSTIK
cat(rep("=", 60), "\n")
## = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = =
par(mfrow = c(2, 2))
# Plot 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)
# Plot 2: Normal Q-Q
qqnorm(residuals_model, pch = 19, main = "Normal Q-Q Plot")
qqline(residuals_model, lwd = 2)
# Plot 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)
# Plot 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))
# ============================================================
# 11. RINGKASAN HASIL ANALISIS
# ============================================================
cat("\n")
cat(rep("=", 60), "\n")
## = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = =
cat("11. RINGKASAN HASIL ANALISIS\n")
## 11. RINGKASAN HASIL ANALISIS
cat(rep("=", 60), "\n")
## = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = =
cat("\nJumlah observasi:", nrow(data), "\n")
##
## Jumlah observasi: 2000
cat("\nModel akhir:", model_compare$Model[idx_final], "\n")
##
## Model akhir: Forward
print(formula(model_final))
## GPA ~ Study_Hours_Per_Day + Extracurricular_Hours_Per_Day + Sleep_Hours_Per_Day
cat("\nKriteria pemilihan model:", kriteria, "\n")
##
## Kriteria pemilihan model: AIC
cat("AIC =", round(AIC(model_final), 4), "\n")
## AIC = -705.9101
cat("BIC =", round(BIC(model_final), 4), "\n")
## BIC = -677.9056
cat("AICc =", round(AICc(model_final), 4), "\n")
## AICc = -705.88
cat("\nUji F: F =", round(f_value, 4), "| p-value =", format(f_pvalue, scientific = TRUE), "\n")
##
## Uji F: F = 783.9318 | p-value = 0e+00
cat("\nR-squared:", round(r_squared, 4), "\n")
##
## R-squared: 0.5409
cat("Adjusted R-squared:", round(adj_r_squared, 4), "\n")
## Adjusted R-squared: 0.5402
cat("\nUji asumsi:\n")
##
## Uji asumsi:
cat("Normalitas (Anderson-Darling) p-value:", format(ad_test$p.value, scientific = TRUE), "\n")
## Normalitas (Anderson-Darling) p-value: 2.711342e-01
cat("Heteroskedastisitas (Breusch-Pagan) p-value:", round(bp_test$p.value, 4), "\n")
## Heteroskedastisitas (Breusch-Pagan) p-value: 0.5816
cat("Autokorelasi (Durbin-Watson) DW:", round(dw_test$statistic, 4), "\n")
## Autokorelasi (Durbin-Watson) DW: 2.0058
cat("VIF maksimum:", vif_maks, "\n")
## VIF maksimum: 1.0008
cat("Linearitas (Ramsey RESET) p-value:", round(reset_test$p.value, 4), "\n")
## Linearitas (Ramsey RESET) p-value: 0.6577
cat("\nDiagnostik:\n")
##
## Diagnostik:
cat("Cook's D maksimum:", round(max(cooks_d), 4), "\n")
## Cook's D maksimum: 0.0115
cat("Leverage maksimum:", round(max(leverage), 4), "\n")
## Leverage maksimum: 0.0047
cat("Jumlah |Studentized Residual| > 3:", sum(abs(std_resid) > 3), "\n")
## Jumlah |Studentized Residual| > 3: 3
cat("\n")
cat(rep("=", 60), "\n")
## = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = =
cat("ANALISIS SELESAI\n")
## ANALISIS SELESAI
cat(rep("=", 60), "\n")
## = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = =