Kelompok 5 :

install.packages(c("car", "corrplot", "nnet"))
library(car)
library(corrplot)
library(nnet)
library(MASS)

Load Dataset

data <- read.csv("mental_health_and_technology_usage_2024.csv")
knitr::kable(head(data))
User_ID Age Gender Technology_Usage_Hours Social_Media_Usage_Hours Gaming_Hours Screen_Time_Hours Mental_Health_Status Stress_Level Sleep_Hours Physical_Activity_Hours Support_Systems_Access Work_Environment_Impact Online_Support_Usage
USER-00001 23 Female 6.57 6.00 0.68 12.36 Good Low 8.01 6.71 No Negative Yes
USER-00002 21 Male 3.01 2.57 3.74 7.61 Poor High 7.28 5.88 Yes Positive No
USER-00003 51 Male 3.04 6.14 1.26 3.16 Fair High 8.04 9.81 No Negative No
USER-00004 25 Female 3.84 4.48 2.59 13.08 Excellent Medium 5.62 5.28 Yes Negative Yes
USER-00005 53 Male 1.20 0.56 0.29 12.63 Good Low 5.55 4.00 No Positive Yes
USER-00006 58 Male 5.59 5.74 0.11 1.34 Poor Low 8.61 6.54 Yes Neutral Yes

Dataset ini membahas hubungan antara penggunaan teknologi sehari-hari (seperti screen time dan media sosial) dengan kondisi kesehatan mental seseorang, seperti tingkat stres, kualitas tidur, dan produktivitas.

Link dataset : https://www.kaggle.com/datasets/waqi786/mental-health-and-technology-usage-dataset

Eksplorasi Data

# Dimensi Data
cat("Dimensi:", nrow(data), "x", ncol(data), "\n")
## Dimensi: 10000 x 14
# Struktur dan Tipe Data
str(data)
## 'data.frame':    10000 obs. of  14 variables:
##  $ User_ID                 : chr  "USER-00001" "USER-00002" "USER-00003" "USER-00004" ...
##  $ Age                     : int  23 21 51 25 53 58 63 51 57 31 ...
##  $ Gender                  : chr  "Female" "Male" "Male" "Female" ...
##  $ Technology_Usage_Hours  : num  6.57 3.01 3.04 3.84 1.2 ...
##  $ Social_Media_Usage_Hours: num  6 2.57 6.14 4.48 0.56 5.74 2.55 4.1 4.11 7.23 ...
##  $ Gaming_Hours            : num  0.68 3.74 1.26 2.59 0.29 0.11 3.79 4.74 0.08 0.81 ...
##  $ Screen_Time_Hours       : num  12.36 7.61 3.16 13.08 12.63 ...
##  $ Mental_Health_Status    : chr  "Good" "Poor" "Fair" "Excellent" ...
##  $ Stress_Level            : chr  "Low" "High" "High" "Medium" ...
##  $ Sleep_Hours             : num  8.01 7.28 8.04 5.62 5.55 8.61 8.61 7.11 7.19 5.09 ...
##  $ Physical_Activity_Hours : num  6.71 5.88 9.81 5.28 4 6.54 1.34 5.27 5.22 0.47 ...
##  $ Support_Systems_Access  : chr  "No" "Yes" "No" "Yes" ...
##  $ Work_Environment_Impact : chr  "Negative" "Positive" "Negative" "Negative" ...
##  $ Online_Support_Usage    : chr  "Yes" "No" "No" "Yes" ...

Interpretasi:

Dataset terdiri dari 10.000 observasi dengan 14 variabel, yang mencakup tipe data campuran yaitu numerik (seperti usia, durasi penggunaan teknologi, waktu tidur, dan aktivitas fisik) serta kategorikal (seperti gender, tingkat stres, dan status kesehatan mental). Selain itu variabel target Mental_Health_Status memiliki lebih dari dua kelas yang bersifat ordinal (Poor < Fair < Good < Excellent), sehingga cocok untuk klasifikasi seperti multinomial logistic regression.

Pre-Processing

1. Ambil 3000 Data

set.seed(123)  

data_sample <- data[sample(nrow(data), 3000), ]

2. Hitung Missing Value

colSums(is.na(data_sample))
##                  User_ID                      Age                   Gender 
##                        0                        0                        0 
##   Technology_Usage_Hours Social_Media_Usage_Hours             Gaming_Hours 
##                        0                        0                        0 
##        Screen_Time_Hours     Mental_Health_Status             Stress_Level 
##                        0                        0                        0 
##              Sleep_Hours  Physical_Activity_Hours   Support_Systems_Access 
##                        0                        0                        0 
##  Work_Environment_Impact     Online_Support_Usage 
##                        0                        0

3. Hapus Data Duplikat

data <- data_sample[!duplicated(data_sample), ]

4. Hapus User_ID

data$User_ID <- NULL

5. Cek Outlier

# Ambil kolom numerik
num_cols <- sapply(data, is.numeric)

# Boxplot semua variabel numerik
boxplot(data[, num_cols], 
        main = "Deteksi Outlier", 
        las = 2)

6. Cek Distribusi Label

table(data$Mental_Health_Status)
## 
## Excellent      Fair      Good      Poor 
##       732       765       748       755

Uji Asumsi

1. Independensi

Catatan: Nilai Durbin-Watson mendekati 2 maka tidak ada autokorelasi, <1.5 maka autokorelasi positif.

model_dummy <- lm(as.numeric(as.factor(Mental_Health_Status)) ~ ., data = data)
# Uji Durbin-Watson
dwt_result <- durbinWatsonTest(model_dummy)
print(dwt_result)
##  lag Autocorrelation D-W Statistic p-value
##    1       0.0101884      1.979004   0.544
##  Alternative hypothesis: rho != 0

Interpretasi: Nilai Durbin Watson Statistic yang diperoleh adalah 1.979 dengan p-value sebesar 0.544. Karena nilai Durbin Watson sangat mendekati angka 2 (berada di dalam rentang aman 1.5 - 2.5) dan p-value > 0.05, kita gagal menolak hipotesis nol. Hal ini berarti tidak ada gejala autokorelasi antar observasi dalam dataset kamu. Asumsi independensi terpenuhi.

2. Non-Multikolinieritas

Catatan: VIF < 10 → aman dan jika VIF < 5 → sangat baik

1. Menggunakan VIF
vif_values <- vif(model_dummy)
print(vif_values)
##                              GVIF Df GVIF^(1/(2*Df))
## Age                      1.005446  1        1.002719
## Gender                   1.009380  2        1.002337
## Technology_Usage_Hours   1.005562  1        1.002777
## Social_Media_Usage_Hours 1.003041  1        1.001519
## Gaming_Hours             1.006673  1        1.003331
## Screen_Time_Hours        1.007137  1        1.003562
## Stress_Level             1.005493  2        1.001370
## Sleep_Hours              1.003109  1        1.001554
## Physical_Activity_Hours  1.005919  1        1.002955
## Support_Systems_Access   1.003110  1        1.001554
## Work_Environment_Impact  1.010694  2        1.002663
## Online_Support_Usage     1.003256  1        1.001626

Interpretasi: Seluruh variabel independen (baik numerik maupun kategorikal) memiliki nilai GVIF di angka sekitar 1.00 hingga 1.01. Karena semua nilai VIF berada jauh di bawah batas maksimal yang diizinkan (VIF < 5 atau VIF < 10) dan plot korelasi juga menunjukkan hubungan antar variabel independen yang sangat lemah, maka terbukti tidak terjadi multikolinieritas. Asumsi ini terpenuhi.

2. Menggunakan Plotting Korelasi
# Mengambil variabel numerik saja
num_data <- data[, sapply(data, is.numeric)]
# Membuat matriks korelasi
cor_matrix <- cor(num_data, use = "complete.obs")
# Plotting korelasi
par(mar = c(5, 4, 10, 2))  

corrplot(cor_matrix, method = "number", type = "upper", 
         tl.cex = 0.6, number.cex = 0.7, 
         tl.srt = 45,
         tl.offset = 1)

title("Matriks Korelasi Fitur Numerik", line = 8)

Interpretasi: Pada visualisasi matriks korelasi, kotak-kotak di luar garis diagonal utama nyaris kosong atau memiliki angka yang sangat mendekati 0, menunjukkan tidak adanya korelasi linear yang kuat antar fitur numerik.

3. Multivariat Outliers

Catatan: Jumlah outlier < 5% data (150) maka data masih dapat digunakan. Dengan menggunakan derajat kebebasan (df) = jumlah variabel numerik, dan alpha = 0.001

# Menghitung Mahalanobis Distance untuk data numerik
m_dist <- mahalanobis(num_data, 
                      center = colMeans(num_data, na.rm = TRUE), 
                      cov = cov(num_data, use = "pairwise.complete.obs"))
# Menentukan nilai kritis (Cut-off) menggunakan distribusi Chi-Square
df <- ncol(num_data)
cutoff <- qchisq(0.999, df)
# Mencari indeks data yang merupakan outlier multivariat
outliers_multi <- which(m_dist > cutoff)
cat("Batas / Cut-off Mahalanobis:", cutoff, "\n")
## Batas / Cut-off Mahalanobis: 24.32189
cat("Jumlah baris yang terdeteksi sebagai outlier multivariat:", length(outliers_multi), "\n")
## Jumlah baris yang terdeteksi sebagai outlier multivariat: 0

Interpretasi: Batas nilai kritis (Cut-off) Mahalanobis yang didapatkan adalah 24.32189. Output menunjukkan bahwa jumlah baris yang terdeteksi melebihi batas cut-off ini adalah 0. Karena tidak ada satupun observasi yang memiliki jarak Mahalanobis di atas cut-off kritis, berarti tidak terdapat data pencilan ekstrem (outlier multivariat) pada sampel 3.000 data yang kamu gunakan. Asumsi ini terpenuhi dengan sempurna.

General Model Ordinal Logistic Regressiion

# Label diurutkan
data$Mental_Health_Status <- ordered(data$Mental_Health_Status, 
                      levels = c("Poor", "Fair", "Good", "Excellent"))

# Model
OLRmodel <- polr(Mental_Health_Status ~ Screen_Time_Hours + Sleep_Hours + 
                 Stress_Level + Technology_Usage_Hours,
                 data = data,
                 Hess = TRUE)

Estimasi Parameter

(OLRestimates <- coef(summary(OLRmodel)))
##                               Value  Std. Error    t value
## Screen_Time_Hours       0.003634888 0.008160177  0.4454422
## Sleep_Hours            -0.003163354 0.022584492 -0.1400675
## Stress_LevelLow         0.043307034 0.080146417  0.5403490
## Stress_LevelMedium      0.111125012 0.080327413  1.3834009
## Technology_Usage_Hours -0.011059656 0.010258443 -1.0781028
## Poor|Fair              -1.102167944 0.183841784 -5.9951983
## Fair|Good               0.015054664 0.182782093  0.0823640
## Good|Excellent          1.120261047 0.184198522  6.0818134

Interpretasi: Karena seluruh nilai p-value dari variabel prediktor bernilai lebih besar dari 0.05, maka tidak ada satupun variabel independen yang berpengaruh secara signifikan terhadap Mental_Health_Status pada dataset sampel yang kamu gunakan.

p <- pnorm(abs(OLRestimates[, "t value"]), lower.tail = FALSE) * 2
p
##      Screen_Time_Hours            Sleep_Hours        Stress_LevelLow 
##           6.560002e-01           8.886066e-01           5.889564e-01 
##     Stress_LevelMedium Technology_Usage_Hours              Poor|Fair 
##           1.665420e-01           2.809879e-01           2.032373e-09 
##              Fair|Good         Good|Excellent 
##           9.343573e-01           1.188308e-09

Interpretasi: Nilai Intercepts (disebut juga \(\alpha\) atau thresholds) menunjukkan ambang batas log-odds yang memisahkan antar kategori pada variabel target yang ordinal saat seluruh variabel prediktor bernilai nol: - Poor|Fair (-1.102): Titik potong (batas) antara kategori Poor dan Fair. Nilai p-value nya sangat kecil (2.03 x 10^{-9}), yang berarti ambang batas ini terdefinisi dengan jelas secara statistik. - Fair|Good (0.015): Titik potong antara kategori Fair dan Good. Nilainya mendekati nol dan p-value nya (0.934) tidak signifikan. - Good|Excellent (1.120): Titik potong antara kategori Good dan Excellent. Nilainya signifikan secara statistik (1.18 x 10^{-9}).

Uji Signifikansi Parameter

1. Uji Serentak

null_model <- polr(Mental_Health_Status ~ 1, data = data, Hess = TRUE)

lrt <- anova(null_model, OLRmodel)
print(lrt)
## Likelihood ratio tests of ordinal regression models
## 
## Response: Mental_Health_Status
##                                                                     Model
## 1                                                                       1
## 2 Screen_Time_Hours + Sleep_Hours + Stress_Level + Technology_Usage_Hours
##   Resid. df Resid. Dev   Test    Df LR stat.  Pr(Chi)
## 1      2997   8316.994                               
## 2      2992   8313.687 1 vs 2     5  3.30685 0.652793
chi_sq <- lrt[2, "LR stat."]
df_lrt <- lrt[2, "   Df"]
p_lrt  <- lrt[2, "Pr(Chi)"]

cat("\n=== Uji Serentak ===\n")
## 
## === Uji Serentak ===
cat("Chi-Square:", chi_sq, "\n")
## Chi-Square: 3.30685
cat("Df        :", df_lrt, "\n")
## Df        : 5
cat("P-value   :", p_lrt, "\n")
## P-value   : 0.652793
if (p_lrt < 0.05) {
  cat("Kesimpulan: Minimal satu variabel berpengaruh signifikan terhadap Mental_Health_Status (tolak H0)\n")
} else {
  cat("Kesimpulan: Tidak ada variabel yang berpengaruh signifikan terhadap Mental_Health_Status (gagal tolak H0)\n")
}
## Kesimpulan: Tidak ada variabel yang berpengaruh signifikan terhadap Mental_Health_Status (gagal tolak H0)

2. Uji Parsial

OLRestimates <- coef(summary(OLRmodel))

p <- pnorm(abs(OLRestimates[, "t value"]), lower.tail = FALSE) * 2

ctable <- cbind(OLRestimates, "p value" = round(p, 4))
print(ctable)
##                               Value  Std. Error    t value p value
## Screen_Time_Hours       0.003634888 0.008160177  0.4454422  0.6560
## Sleep_Hours            -0.003163354 0.022584492 -0.1400675  0.8886
## Stress_LevelLow         0.043307034 0.080146417  0.5403490  0.5890
## Stress_LevelMedium      0.111125012 0.080327413  1.3834009  0.1665
## Technology_Usage_Hours -0.011059656 0.010258443 -1.0781028  0.2810
## Poor|Fair              -1.102167944 0.183841784 -5.9951983  0.0000
## Fair|Good               0.015054664 0.182782093  0.0823640  0.9344
## Good|Excellent          1.120261047 0.184198522  6.0818134  0.0000
cat("Interpretasi Uji Parsial")
## Interpretasi Uji Parsial
predictor_rows <- rownames(ctable)[!grepl("\\|", rownames(ctable))]

for (var in predictor_rows) {
  pval <- ctable[var, "p value"]
  sig <- ifelse(pval < 0.05, "SIGNIFIKAN ✓", "tidak signifikan ✗")
  cat(sprintf("%-30s p-value = %.4f → %s\n", var, pval, sig))
}
## Screen_Time_Hours              p-value = 0.6560 → tidak signifikan ✗
## Sleep_Hours                    p-value = 0.8886 → tidak signifikan ✗
## Stress_LevelLow                p-value = 0.5890 → tidak signifikan ✗
## Stress_LevelMedium             p-value = 0.1665 → tidak signifikan ✗
## Technology_Usage_Hours         p-value = 0.2810 → tidak signifikan ✗

3. Odds Ratio

OR <- exp(coef(OLRmodel))

# Confidence Interval 95%
ci <- exp(confint(OLRmodel))
## Waiting for profiling to be done...
or_table <- cbind("Odds Ratio" = round(OR, 4),
                  "CI 2.5%"    = round(ci[, 1], 4),
                  "CI 97.5%"   = round(ci[, 2], 4))
print(or_table)
##                        Odds Ratio CI 2.5% CI 97.5%
## Screen_Time_Hours          1.0036  0.9877   1.0198
## Sleep_Hours                0.9968  0.9537   1.0420
## Stress_LevelLow            1.0443  0.8925   1.2219
## Stress_LevelMedium         1.1175  0.9548   1.3082
## Technology_Usage_Hours     0.9890  0.9693   1.0091
cat("\n=== Interpretasi Odds Ratio ===\n")
## 
## === Interpretasi Odds Ratio ===
for (var in names(OR)) {
  arah <- ifelse(OR[var] > 1, "meningkatkan", "menurunkan")
  cat(sprintf("%-30s OR = %.4f → %s peluang Mental Health lebih tinggi\n", 
              var, OR[var], arah))
}
## Screen_Time_Hours              OR = 1.0036 → meningkatkan peluang Mental Health lebih tinggi
## Sleep_Hours                    OR = 0.9968 → menurunkan peluang Mental Health lebih tinggi
## Stress_LevelLow                OR = 1.0443 → meningkatkan peluang Mental Health lebih tinggi
## Stress_LevelMedium             OR = 1.1175 → meningkatkan peluang Mental Health lebih tinggi
## Technology_Usage_Hours         OR = 0.9890 → menurunkan peluang Mental Health lebih tinggi

4. Goddess Of Fit

# McFadden R2
null_loglik  <- logLik(null_model)
full_loglik  <- logLik(OLRmodel)

mcfadden_r2 <- as.numeric(1 - (full_loglik / null_loglik))

# Cox & Snell R2
n <- nrow(data)
cox_snell_r2 <- as.numeric(1 - exp((2/n) * (null_loglik - full_loglik)))

# Nagelkerke R2
nagelkerke_r2 <- cox_snell_r2 / (1 - exp((2/n) * null_loglik))

cat("Pseudo R-Squared")
## Pseudo R-Squared
cat(sprintf("McFadden R²    : %.4f\n", mcfadden_r2))
## McFadden R²    : 0.0004
cat(sprintf("Cox & Snell R² : %.4f\n", cox_snell_r2))
## Cox & Snell R² : 0.0011
cat(sprintf("Nagelkerke R²  : %.4f\n", nagelkerke_r2))
## Nagelkerke R²  : 0.0012