Analisis MANCOVA
Kelompok 7
2025-06-05
Dosen Pembimbing: Ike Fitriyaningsih, M.Si
NIDN : 0109049001
Kelas : 2023E
Pendahuluan
Dokumen ini memuat hasil analisis MANCOVA pada dataset Adult
Income dari UCI. Tujuan dari analisis ini adalah untuk mengetahui
pengaruh variabel age, occupation, dan
income terhadap dua variabel dependen yaitu
capital.gain dan capital.loss.
Seluruh proses analisis mencakup: pembersihan data, pengujian asumsi MANCOVA (normalitas, homogenitas varians, dll), transformasi data, dan interpretasi hasil.
Load library yang diperlukan
library(tidyverse)## Warning: package 'tidyverse' was built under R version 4.4.3## Warning: package 'ggplot2' was built under R version 4.4.3## Warning: package 'tibble' was built under R version 4.4.3## Warning: package 'tidyr' was built under R version 4.4.3## Warning: package 'readr' was built under R version 4.4.3## Warning: package 'purrr' was built under R version 4.4.3## Warning: package 'dplyr' was built under R version 4.4.3## Warning: package 'stringr' was built under R version 4.4.3## Warning: package 'forcats' was built under R version 4.4.3## Warning: package 'lubridate' was built under R version 4.4.3## ── Attaching core tidyverse packages ──────────────────────── tidyverse 2.0.0 ──
## ✔ dplyr 1.1.4 ✔ readr 2.1.5
## ✔ forcats 1.0.0 ✔ stringr 1.5.1
## ✔ ggplot2 3.5.1 ✔ tibble 3.2.1
## ✔ lubridate 1.9.4 ✔ tidyr 1.3.1
## ✔ purrr 1.0.4
## ── Conflicts ────────────────────────────────────────── tidyverse_conflicts() ──
## ✖ dplyr::filter() masks stats::filter()
## ✖ dplyr::lag() masks stats::lag()
## ℹ Use the conflicted package (<http://conflicted.r-lib.org/>) to force all conflicts to become errorslibrary(car)## Warning: package 'car' was built under R version 4.4.3## Loading required package: carData## Warning: package 'carData' was built under R version 4.4.3##
## Attaching package: 'car'
##
## The following object is masked from 'package:dplyr':
##
## recode
##
## The following object is masked from 'package:purrr':
##
## somelibrary(MVN)## Warning: package 'MVN' was built under R version 4.4.3library(biotools)## Warning: package 'biotools' was built under R version 4.4.3## Loading required package: MASS
##
## Attaching package: 'MASS'
##
## The following object is masked from 'package:dplyr':
##
## select
##
## ---
## biotools version 4.3library(ggpubr)## Warning: package 'ggpubr' was built under R version 4.4.3Memanggil pustaka-pustaka yang dibutuhkan seperti tidyverse (manipulasi data), car (analisis ANOVA/MANCOVA), MVN (uji normalitas multivariat), biotools (uji Box’s M), dan ggpubr (grafik).
Load Data
df <- read.csv("C:/Users/HP 430 G5/OneDrive/Dokumen/Semester 4/Analisis Multivariat/adult.csv", stringsAsFactors = TRUE)Cek struktur data dan Missing value
cat("Struktur data awal:\n")## Struktur data awal:str(df)## 'data.frame': 32561 obs. of 15 variables:
## $ age : int 90 82 66 54 41 34 38 74 68 41 ...
## $ workclass : Factor w/ 9 levels "?","Federal-gov",..: 1 5 1 5 5 5 5 8 2 5 ...
## $ fnlwgt : int 77053 132870 186061 140359 264663 216864 150601 88638 422013 70037 ...
## $ education : Factor w/ 16 levels "10th","11th",..: 12 12 16 6 16 12 1 11 12 16 ...
## $ education.num : int 9 9 10 4 10 9 6 16 9 10 ...
## $ marital.status: Factor w/ 7 levels "Divorced","Married-AF-spouse",..: 7 7 7 1 6 1 6 5 1 5 ...
## $ occupation : Factor w/ 15 levels "?","Adm-clerical",..: 1 5 1 8 11 9 2 11 11 4 ...
## $ relationship : Factor w/ 6 levels "Husband","Not-in-family",..: 2 2 5 5 4 5 5 3 2 5 ...
## $ race : Factor w/ 5 levels "Amer-Indian-Eskimo",..: 5 5 3 5 5 5 5 5 5 5 ...
## $ sex : Factor w/ 2 levels "Female","Male": 1 1 1 1 1 1 2 1 1 2 ...
## $ capital.gain : int 0 0 0 0 0 0 0 0 0 0 ...
## $ capital.loss : int 4356 4356 4356 3900 3900 3770 3770 3683 3683 3004 ...
## $ hours.per.week: int 40 18 40 40 40 45 40 20 40 60 ...
## $ native.country: Factor w/ 42 levels "?","Cambodia",..: 40 40 40 40 40 40 40 40 40 1 ...
## $ income : Factor w/ 2 levels "<=50K",">50K": 1 1 1 1 1 1 1 2 1 2 ...cat("\n5 baris pertama dari data:\n")##
## 5 baris pertama dari data:head(df, 5)## age workclass fnlwgt education education.num marital.status
## 1 90 ? 77053 HS-grad 9 Widowed
## 2 82 Private 132870 HS-grad 9 Widowed
## 3 66 ? 186061 Some-college 10 Widowed
## 4 54 Private 140359 7th-8th 4 Divorced
## 5 41 Private 264663 Some-college 10 Separated
## occupation relationship race sex capital.gain capital.loss
## 1 ? Not-in-family White Female 0 4356
## 2 Exec-managerial Not-in-family White Female 0 4356
## 3 ? Unmarried Black Female 0 4356
## 4 Machine-op-inspct Unmarried White Female 0 3900
## 5 Prof-specialty Own-child White Female 0 3900
## hours.per.week native.country income
## 1 40 United-States <=50K
## 2 18 United-States <=50K
## 3 40 United-States <=50K
## 4 40 United-States <=50K
## 5 40 United-States <=50Kcat("\nJumlah missing value per kolom:\n")##
## Jumlah missing value per kolom:colSums(is.na(df) | df == " ?")## age workclass fnlwgt education education.num
## 0 0 0 0 0
## marital.status occupation relationship race sex
## 0 0 0 0 0
## capital.gain capital.loss hours.per.week native.country income
## 0 0 0 0 0Tangani missing value
df[df == " ?"] <- NA
cat("\nJumlah missing value setelah konversi ' ?' ke NA:\n")##
## Jumlah missing value setelah konversi ' ?' ke NA:colSums(is.na(df))## age workclass fnlwgt education education.num
## 0 0 0 0 0
## marital.status occupation relationship race sex
## 0 0 0 0 0
## capital.gain capital.loss hours.per.week native.country income
## 0 0 0 0 0Mengganti tanda ” ?” menjadi NA agar bisa dibersihkan dan dianalisis secara statistik.
Sampling dan hapus NA
df <- df[1:4000, ]
df <- na.omit(df)
cat("\nDimensi data setelah sampling dan menghapus NA:\n")##
## Dimensi data setelah sampling dan menghapus NA:dim(df)## [1] 4000 15mengambil 4000 baris pertama dari dataset untuk analisis dan menghapus baris yang memiliki nilai NA (missing value) dengan fungsi na.omit(). mengambil sampel 4000 data pertama, kita mendapatkan ukuran data yang cukup representatif untuk analisis, tapi tetap efisien dalam hal waktu dan sumber daya komputasi.
Pilih variabel untuk MANCOVA
data_mancova <- df[, c("age", "hours.per.week", "occupation", "income", "capital.gain", "capital.loss")]
cat("\nStruktur data_mancova:\n")##
## Struktur data_mancova:str(data_mancova)## 'data.frame': 4000 obs. of 6 variables:
## $ age : int 90 82 66 54 41 34 38 74 68 41 ...
## $ hours.per.week: int 40 18 40 40 40 45 40 20 40 60 ...
## $ occupation : Factor w/ 15 levels "?","Adm-clerical",..: 1 5 1 8 11 9 2 11 11 4 ...
## $ income : Factor w/ 2 levels "<=50K",">50K": 1 1 1 1 1 1 1 2 1 2 ...
## $ capital.gain : int 0 0 0 0 0 0 0 0 0 0 ...
## $ capital.loss : int 4356 4356 4356 3900 3900 3770 3770 3683 3683 3004 ...capital.gain dan capital.loss → Ini adalah variabel dependen dalam analisis MANCOVA. Kita ingin melihat apakah dan bagaimana dua jenis modal finansial ini dipengaruhi oleh variabel lain secara simultan.
age → Dimasukkan sebagai kovariat, karena umur bisa memengaruhi jumlah capital gain/loss yang diperoleh seseorang, tapi bukan faktor eksperimen utama.
occupation dan income → Dimasukkan sebagai variabel independen (faktor), karena kita ingin mengetahui apakah pekerjaan dan tingkat pendapatan seseorang berpengaruh terhadap capital gain dan loss mereka.
Konversi faktor dan ringkasan
data_mancova$occupation <- factor(data_mancova$occupation)
data_mancova$income <- factor(data_mancova$income)
cat("\nRingkasan statistik awal:\n")##
## Ringkasan statistik awal:summary(data_mancova)## age hours.per.week occupation income
## Min. :17.00 Min. : 1.0 Exec-managerial :826 <=50K:1550
## 1st Qu.:34.00 1st Qu.:40.0 Prof-specialty :823 >50K :2450
## Median :42.00 Median :40.0 Sales :493
## Mean :43.31 Mean :43.9 Craft-repair :489
## 3rd Qu.:51.00 3rd Qu.:50.0 Adm-clerical :329
## Max. :90.00 Max. :99.0 Machine-op-inspct:182
## (Other) :858
## capital.gain capital.loss
## Min. : 0 Min. : 0.0
## 1st Qu.: 0 1st Qu.: 0.0
## Median : 3464 Median : 0.0
## Mean : 8688 Mean : 710.7
## 3rd Qu.: 7688 3rd Qu.:1762.0
## Max. :99999 Max. :4356.0
## Memastikan struktur data sudah sesuai (tipe data dan distribusi kategori benar). Memberi gambaran awal mengenai penyebaran data yang akan dianalisis lebih lanjut.
TAHAP ANALISIS
Uji normalitas Multivariat sebelum transformasi
normal_pre <- mvn(data_mancova[, c("capital.gain", "capital.loss")], mvnTest = "hz")
normal_pre$multivariateNormality## Test HZ p value MVN
## 1 Henze-Zirkler 976.6442 0 NOgunakan fungsi mvn() dari package MVN dengan metode Henze-Zirkler, yang menguji apakah gabungan dari variabel capital.gain dan capital.loss mengikuti distribusi normal multivariat.
Nilai HZ tinggi dan p-value = 0, menunjukkan penolakan terhadap normalitas multivariat. Artinya, data capital.gain dan capital.loss secara gabungan tidak terdistribusi normal.
Lakukan Transformasi log dan uji1p lagi
data_mancova$capital.gain <- log1p(data_mancova$capital.gain)
data_mancova$capital.loss <- log1p(data_mancova$capital.loss)
normal_post <- mvn(data_mancova[, c("capital.gain", "capital.loss")], mvnTest = "hz")
normal_post$multivariateNormality## Test HZ p value MVN
## 1 Henze-Zirkler 562.4332 0 NOKarena data tidak normal, kita lakukan transformasi log1p pada kedua variabel Meskipun nilai HZ menurun setelah transformasi, p-value tetap 0 dan hasilnya masih tidak normal. Transformasi logaritmik membantu mengurangi skewness, tapi belum cukup untuk memenuhi asumsi normalitas multivariat.
Uji Homogenitas Kovarians menggunkan Box-M
data_mancova$group <- interaction(data_mancova$occupation, data_mancova$income)
group_count <- table(data_mancova$group)
valid_groups <- names(group_count[group_count >= 3])
data_valid <- data_mancova[data_mancova$group %in% valid_groups, ]
boxm_result <- boxM(data_valid[, c("capital.gain", "capital.loss")], data_valid$group)
boxm_result##
## Box's M-test for Homogeneity of Covariance Matrices
##
## data: data_valid[, c("capital.gain", "capital.loss")]
## Chi-Sq (approx.) = 1534.4, df = 78, p-value < 2.2e-16Uji ini mengevaluasi apakah matriks kovarians antar grup sama atau homogen Nilai p-value jauh lebih kecil dari 0.05, maka kita menolak H0 (hipotesis nol bahwa semua matriks kovarians antar grup adalah sama). Artinya, terdapat perbedaan signifikan pada matriks kovarians antar kombinasi occupation dan income.
Uji Homogenitas Varians dengan Uji Levene
leveneTest(capital.gain ~ occupation * income, data = data_mancova)## Levene's Test for Homogeneity of Variance (center = median)
## Df F value Pr(>F)
## group 28 1.6961 0.01254 *
## 3971
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1leveneTest(capital.loss ~ occupation * income, data = data_mancova)## Levene's Test for Homogeneity of Variance (center = median)
## Df F value Pr(>F)
## group 28 3.7766 8.638e-11 ***
## 3971
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1Menguji apakah varians antar kelompok (gabungan occupation * income) untuk masing-masing variabel dependen (capital.gain dan capital.loss) adalah homogen.
Untuk capital.gain: F value = 1.6961 p-value = 0.01254 → Signifikan
Untuk capital.loss: F value = 3.7766 p-value = 8.638e-11 → Sangat signifikan
Kedua p-value < 0.05 → Asumsi homogenitas varians TIDAK terpenuhi untuk kedua variabel dependen. Artinya, terdapat perbedaan varians antara kelompok occupation * income pada capital.gain dan capital.loss.
Uji linearitas dengan scatter plot
scatterplotMatrix(data_mancova[, c("age", "capital.gain", "capital.loss")],
smooth = FALSE)
age vs capital.gain Tampak pola menyebar horizontal dengan garis tren
positif tapi tidak kuat. Kesimpulan: Ada indikasi hubungan linier lemah,
namun kemungkinan besar data ini tidak sepenuhnya linier.
age vs capital.loss Pola sebaran cukup horizontal, garis tren cenderung menurun. Kesimpulan: Ada sedikit indikasi hubungan negatif, namun juga tidak sepenuhnya linier.
capital.gain vs capital.loss
Sangat terlihat berpola L terbalik, menunjukkan non-linear relationship. Garis tren menurun tajam, menandakan bahwa jika salah satu tinggi, yang lain cenderung nol. Ini logis karena data sering menunjukkan orang hanya mendapat gain atau loss, jarang dua-duanya.
Hubungan antara variabel kovariat (age) dengan variabel dependen (capital.gain, capital.loss) tidak sepenuhnya linier. Pola sebar dan fitting menunjukkan bahwa asumsi linearitas tidak terpenuhi secara kuat.
Uji Multikolinearitas
lm_check <- lm(age ~ occupation + income, data = data_mancova)
vif(lm_check)## GVIF Df GVIF^(1/(2*Df))
## occupation 1.160045 14 1.005316
## income 1.160045 1 1.077054Hasil uji multikolinearitas menunjukkan bahwa nilai GVIF^(1/(2×Df)) untuk variabel occupation adalah 1.005 dan untuk income adalah 1.077, yang keduanya jauh di bawah ambang batas 2. Ini mengindikasikan bahwa tidak terdapat masalah multikolinearitas antara variabel-variabel independen dalam model, sehingga aman digunakan untuk analisis lebih lanjut
Model MANCOVA
model <- manova(cbind(capital.gain, capital.loss) ~ age + occupation + income, data = data_mancova)
summary(model, test = "Pillai")## Df Pillai approx F num Df den Df Pr(>F)
## age 1 0.030824 63.32 2 3982 < 2.2e-16 ***
## occupation 14 0.095792 14.31 28 7966 < 2.2e-16 ***
## income 1 0.188169 461.48 2 3982 < 2.2e-16 ***
## Residuals 3983
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1Semua variabel (age, occupation, income) berpengaruh signifikan terhadap gabungan variabel dependen (trans_hours, trans_edu). p-value < 0.001 pada semuanya berarti pengaruhnya sangat signifikan secara statistik. Nilai Pillai’s trace mengukur kekuatan pengaruh: - income punya pengaruh terbesar (Pillai: 0.188) - disusul occupation dan age.
Pillai’s Trace digunakan dalam MANCOVA karena merupakan statistik yang paling robust terhadap pelanggaran asumsi, seperti normalitas multivariat dan homogenitas kovarian. Dibandingkan metode lain (seperti Wilks’ Lambda), Pillai lebih stabil saat data tidak seimbang dan lebih konservatif, sehingga mengurangi risiko kesalahan tipe I. Oleh karena itu, dalam kondisi data yang kompleks dan tidak sepenuhnya memenuhi asumsi klasik, Pillai menjadi pilihan yang paling aman untuk interpretasi hasil.
summary.aov(model)## Response capital.gain :
## Df Sum Sq Mean Sq F value Pr(>F)
## age 1 1025 1024.7 55.8846 9.397e-14 ***
## occupation 14 552 39.4 2.1493 0.00758 **
## income 1 3292 3292.3 179.5522 < 2.2e-16 ***
## Residuals 3983 73034 18.3
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Response capital.loss :
## Df Sum Sq Mean Sq F value Pr(>F)
## age 1 487 487.30 37.6833 9.134e-10 ***
## occupation 14 140 10.02 0.7748 0.6979
## income 1 1092 1091.54 84.4104 < 2.2e-16 ***
## Residuals 3983 51506 12.93
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1Hipotesis Nol (H₀): Setelah mengontrol pengaruh age, tidak terdapat perbedaan yang signifikan pada rata-rata gabungan antara capital.loss dan capital.gain berdasarkan kombinasi dari occupation, dan income. Hipotesis Alternatif (H₁): Setelah mengontrol pengaruh age, terdapat perbedaan yang signifikan pada rata-rata gabungan antara capital.loss dan capital.gain berdasarkan paling tidak satu kelompok dari occupation, atau income.
Semua variabel bermakna secara multivariat, dengan nilai p < 0.001. Maka, tolak H₀. Setelah mengontrol age, terdapat perbedaan yang signifikan dalam rata-rata gabungan antara capital.gain dan capital.loss berdasarkan paling tidak satu dari occupation, atau income.