Modul 2 Implementasi Multivariate Analysis of Variance (MANOVA) dan Multivariate Analysis of Covariance (MANCOVA) pada Parameter Fungsi Hati Berdasarkan HCV Dataset.
Ananta Putri Clara (24031554031), Ayu Anggraini (24031554054), Faizah Khasifa Salsabila (24031554056)
LOAD LIBRARY
install.packages(c("tidyverse", "car"))## Installing packages into '/cloud/lib/x86_64-pc-linux-gnu-library/4.5'
## (as 'lib' is unspecified)Pada tahap awal, dilakukan instalasi package tidyverse dan car agar seluruh kebutuhan analisis, mulai dari pengolahan data hingga pengujian statistik, dapat terpenuhi.
library(tidyverse)## ── Attaching core tidyverse packages ──────────────────────── tidyverse 2.0.0 ──
## ✔ dplyr 1.2.0 ✔ readr 2.2.0
## ✔ forcats 1.0.1 ✔ stringr 1.6.0
## ✔ ggplot2 4.0.2 ✔ tibble 3.3.1
## ✔ lubridate 1.9.5 ✔ tidyr 1.3.2
## ✔ purrr 1.2.1
## ── 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 errorsSetelah itu, tidyverse dipanggil ke dalam R sehingga fungsi-fungsi untuk manipulasi data dan visualisasi dapat langsung digunakan.
library(car)## Loading required package: carData##
## Attaching package: 'car'## The following object is masked from 'package:dplyr':
##
## recode## The following object is masked from 'package:purrr':
##
## someSelanjutnya, package car juga diaktifkan untuk mendukung pengujian statistik, khususnya dalam uji asumsi seperti homogenitas varians.
LOAD DATA
data <- read.csv("hcvdata.csv")Berikutnya, dataset dibaca dari file hcvdata.csv dan dimasukkan ke dalam variabel data untuk diproses lebih lanjut.
data <- data[, -1]Setelah data dimuat, kolom pertama dihapus karena umumnya hanya berisi indeks yang tidak diperlukan dalam analisis.
head(data)## Category Age Sex ALB ALP ALT AST BIL CHE CHOL CREA GGT PROT
## 1 0=Blood Donor 32 m 38.5 52.5 7.7 22.1 7.5 6.93 3.23 106 12.1 69.0
## 2 0=Blood Donor 32 m 38.5 70.3 18.0 24.7 3.9 11.17 4.80 74 15.6 76.5
## 3 0=Blood Donor 32 m 46.9 74.7 36.2 52.6 6.1 8.84 5.20 86 33.2 79.3
## 4 0=Blood Donor 32 m 43.2 52.0 30.6 22.6 18.9 7.33 4.74 80 33.8 75.7
## 5 0=Blood Donor 32 m 39.2 74.1 32.6 24.8 9.6 9.15 4.32 76 29.9 68.7
## 6 0=Blood Donor 32 m 41.6 43.3 18.5 19.7 12.3 9.92 6.05 111 91.0 74.0Untuk memastikan data telah terbaca dengan benar, ditampilkan beberapa baris awal sebagai gambaran isi dataset.
str(data)## 'data.frame': 615 obs. of 13 variables:
## $ Category: chr "0=Blood Donor" "0=Blood Donor" "0=Blood Donor" "0=Blood Donor" ...
## $ Age : int 32 32 32 32 32 32 32 32 32 32 ...
## $ Sex : chr "m" "m" "m" "m" ...
## $ ALB : num 38.5 38.5 46.9 43.2 39.2 41.6 46.3 42.2 50.9 42.4 ...
## $ ALP : num 52.5 70.3 74.7 52 74.1 43.3 41.3 41.9 65.5 86.3 ...
## $ ALT : num 7.7 18 36.2 30.6 32.6 18.5 17.5 35.8 23.2 20.3 ...
## $ AST : num 22.1 24.7 52.6 22.6 24.8 19.7 17.8 31.1 21.2 20 ...
## $ BIL : num 7.5 3.9 6.1 18.9 9.6 12.3 8.5 16.1 6.9 35.2 ...
## $ CHE : num 6.93 11.17 8.84 7.33 9.15 ...
## $ CHOL : num 3.23 4.8 5.2 4.74 4.32 6.05 4.79 4.6 4.1 4.45 ...
## $ CREA : num 106 74 86 80 76 111 70 109 83 81 ...
## $ GGT : num 12.1 15.6 33.2 33.8 29.9 91 16.9 21.5 13.7 15.9 ...
## $ PROT : num 69 76.5 79.3 75.7 68.7 74 74.5 67.1 71.3 69.9 ...Kemudian, struktur data ditinjau untuk mengetahui tipe setiap variabel serta bentuk keseluruhan dataset.
PREPROCESSING
data$Sex <- as.factor(data$Sex)
data$Category <- as.factor(data$Category)Selanjutnya, variabel Sex diubah menjadi tipe faktor agar sesuai sebagai data kategorikal dalam analisis lalu pada variabel Category sehingga dapat diperlakukan sebagai variabel kategori.
colSums(is.na(data))## Category Age Sex ALB ALP ALT AST BIL
## 0 0 0 1 18 1 0 0
## CHE CHOL CREA GGT PROT
## 0 10 0 0 1Setelah itu, dilakukan pengecekan terhadap missing value untuk melihat apakah terdapat data yang kosong pada setiap kolom.
data_bersih <- na.omit(data)Jika ditemukan nilai yang hilang, data tersebut kemudian dibersihkan dengan menghapus baris yang mengandung missing value sehingga diperoleh dataset yang siap dianalisis.
STATISTIK DESKRIPTIF
summary(data_bersih[, c("Age","ALB","ALP","ALT","AST","BIL",
"CHE","CHOL","CREA","GGT","PROT")])## Age ALB ALP ALT
## Min. :23.00 Min. :14.90 Min. : 11.30 Min. : 0.90
## 1st Qu.:39.00 1st Qu.:38.80 1st Qu.: 52.50 1st Qu.: 16.40
## Median :47.00 Median :41.90 Median : 66.20 Median : 22.70
## Mean :47.42 Mean :41.62 Mean : 68.12 Mean : 26.58
## 3rd Qu.:54.00 3rd Qu.:45.10 3rd Qu.: 79.90 3rd Qu.: 31.90
## Max. :77.00 Max. :82.20 Max. :416.60 Max. :325.30
## AST BIL CHE CHOL
## Min. : 10.60 Min. : 0.80 Min. : 1.420 Min. :1.430
## 1st Qu.: 21.50 1st Qu.: 5.20 1st Qu.: 6.930 1st Qu.:4.620
## Median : 25.70 Median : 7.10 Median : 8.260 Median :5.310
## Mean : 33.77 Mean : 11.02 Mean : 8.204 Mean :5.391
## 3rd Qu.: 31.70 3rd Qu.: 11.00 3rd Qu.: 9.570 3rd Qu.:6.080
## Max. :324.00 Max. :209.00 Max. :16.410 Max. :9.670
## CREA GGT PROT
## Min. : 8.00 Min. : 4.5 Min. :44.80
## 1st Qu.: 68.00 1st Qu.: 15.6 1st Qu.:69.30
## Median : 77.00 Median : 22.8 Median :72.10
## Mean : 81.67 Mean : 38.2 Mean :71.89
## 3rd Qu.: 89.00 3rd Qu.: 37.6 3rd Qu.:75.20
## Max. :1079.10 Max. :650.9 Max. :86.50table(data_bersih$Sex)##
## f m
## 226 363table(data_bersih$Category)##
## 0=Blood Donor 0s=suspect Blood Donor 1=Hepatitis
## 526 7 20
## 2=Fibrosis 3=Cirrhosis
## 12 24Selanjutnya, ditampilkan ringkasan statistik dari variabel numerik guna memperoleh gambaran awal mengenai distribusi dan karakteristik data, kemudian distribusi data berdasarkan jenis kelamin ditinjau untuk mengetahui jumlah masing-masing kategori, lalu hal yang sama juga dilakukan pada variabel kategori pasien untuk melihat sebaran data pada tiap kelompok.
VISUALISASI
data_bersih %>%
select(ALB, ALP, ALT, AST, BIL, CHE, CHOL, CREA, GGT, PROT) %>%
pivot_longer(cols = everything(), names_to = "Variabel", values_to = "Nilai") %>%
ggplot(aes(x = Nilai)) +
geom_histogram(bins = 20, fill = "steelblue", color = "black") +
facet_wrap(~Variabel, scales = "free") +
labs(title = "Distribusi Parameter Fungsi Hati")
Berikutnya, dilakukan visualisasi distribusi data dengan mengubah
dataset ke format panjang terlebih dahulu, kemudian ditampilkan dalam
bentuk histogram untuk masing-masing parameter fungsi hati sehingga pola
distribusi dapat diamati dengan lebih jelas.
data_bersih %>%
select(Category, ALB, ALP, ALT, AST, BIL, CHE, CHOL, CREA, GGT, PROT) %>%
pivot_longer(cols = -Category, names_to = "Variabel", values_to = "Nilai") %>%
ggplot(aes(x = Category, y = Nilai, fill = Category)) +
geom_boxplot() +
facet_wrap(~Variabel, scales = "free") +
labs(title = "Boxplot Parameter Fungsi Hati berdasarkan Kategori Pasien")
Selanjutnya, dibuat boxplot untuk membandingkan distribusi setiap
parameter fungsi hati berdasarkan kategori pasien sehingga perbedaan
antar kelompok dapat terlihat dengan lebih jelas.
UJI ASUMSI
shapiro.test(data_bersih$ALB)##
## Shapiro-Wilk normality test
##
## data: data_bersih$ALB
## W = 0.92906, p-value = 4.502e-16shapiro.test(data_bersih$AST)##
## Shapiro-Wilk normality test
##
## data: data_bersih$AST
## W = 0.43504, p-value < 2.2e-16Setelah visualisasi, dilakukan pengujian normalitas pada variabel ALB untuk melihat apakah data mengikuti distribusi normal, dan kemudian pengujian yang sama juga dilakukan pada variabel AST sebagai representasi dari variabel lainnya.
leveneTest(ALB ~ Category, data = data_bersih)## Levene's Test for Homogeneity of Variance (center = median)
## Df F value Pr(>F)
## group 4 1.2525 0.2876
## 584leveneTest(AST ~ Category, data = data_bersih)## Levene's Test for Homogeneity of Variance (center = median)
## Df F value Pr(>F)
## group 4 68.273 < 2.2e-16 ***
## 584
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1Kemudian, diuji apakah varians data pada variabel ALB seragam antar kategori sebagai bagian dari uji homogenitas, langkah ini diulang pada variabel AST untuk memastikan asumsi homogenitas terpenuhi.
MANOVA
manova_model <- manova(
cbind(ALB, ALP, ALT, AST, BIL, CHE, CHOL, CREA, GGT, PROT) ~ Category + Sex,
data = data_bersih
)Selanjutnya, dibangun model MANOVA dengan beberapa parameter fungsi hati sebagai variabel dependen dan kategori serta jenis kelamin sebagai variabel independen.
summary(manova_model, test = "Pillai")## Df Pillai approx F num Df den Df Pr(>F)
## Category 4 1.38633 30.605 40 2308 < 2.2e-16 ***
## Sex 1 0.15728 10.712 10 574 < 2.2e-16 ***
## Residuals 583
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1summary(manova_model, test = "Wilks")## Df Wilks approx F num Df den Df Pr(>F)
## Category 4 0.12220 40.347 40 2178.4 < 2.2e-16 ***
## Sex 1 0.84272 10.712 10 574.0 < 2.2e-16 ***
## Residuals 583
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1summary.aov(manova_model)## Response ALB :
## Df Sum Sq Mean Sq F value Pr(>F)
## Category 4 4587.0 1146.75 46.594 < 2e-16 ***
## Sex 1 585.1 585.08 23.773 1.4e-06 ***
## Residuals 583 14348.5 24.61
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Response ALP :
## Df Sum Sq Mean Sq F value Pr(>F)
## Category 4 51143 12785.7 21.6735 <2e-16 ***
## Sex 1 9 9.4 0.0159 0.8996
## Residuals 583 343926 589.9
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Response ALT :
## Df Sum Sq Mean Sq F value Pr(>F)
## Category 4 49478 12369.4 36.147 < 2.2e-16 ***
## Sex 1 6962 6962.2 20.346 7.819e-06 ***
## Residuals 583 199499 342.2
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Response AST :
## Df Sum Sq Mean Sq F value Pr(>F)
## Category 4 276232 69058 113.6325 < 2.2e-16 ***
## Sex 1 4638 4638 7.6311 0.005918 **
## Residuals 583 354307 608
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Response BIL :
## Df Sum Sq Mean Sq F value Pr(>F)
## Category 4 62883 15720.8 80.5338 < 2.2e-16 ***
## Sex 1 1468 1468.0 7.5202 0.006289 **
## Residuals 583 113806 195.2
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Response CHE :
## Df Sum Sq Mean Sq F value Pr(>F)
## Category 4 590.54 147.636 40.391 < 2.2e-16 ***
## Sex 1 101.35 101.346 27.727 1.968e-07 ***
## Residuals 583 2130.98 3.655
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Response CHOL :
## Df Sum Sq Mean Sq F value Pr(>F)
## Category 4 69.41 17.3531 14.8805 1.381e-11 ***
## Sex 1 0.14 0.1431 0.1227 0.7262
## Residuals 583 679.87 1.1662
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Response CREA :
## Df Sum Sq Mean Sq F value Pr(>F)
## Category 4 138855 34714 15.154 8.547e-12 ***
## Sex 1 36917 36917 16.116 6.737e-05 ***
## Residuals 583 1335497 2291
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Response GGT :
## Df Sum Sq Mean Sq F value Pr(>F)
## Category 4 455756 113939 52.5243 < 2e-16 ***
## Sex 1 13428 13428 6.1902 0.01312 *
## Residuals 583 1264681 2169
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Response PROT :
## Df Sum Sq Mean Sq F value Pr(>F)
## Category 4 2705.6 676.41 28.0513 <2e-16 ***
## Sex 1 59.4 59.35 2.4614 0.1172
## Residuals 583 14058.0 24.11
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1Hasil analisis kemudian ditinjau menggunakan statistik Pillai Trace untuk melihat pengaruh variabel independen secara simultan, selain itu, digunakan juga Wilks Lambda sebagai pembanding dalam menguji signifikansi model, untuk memperdalam analisis, dilihat pula hasil uji pada masing-masing variabel secara terpisah.
MANCOVA
mancova_model <- manova(
cbind(ALB, ALP, ALT, AST, BIL, CHE, CHOL, CREA, GGT, PROT) ~ Category + Sex + Age,
data = data_bersih
)Terakhir, model dikembangkan menjadi MANCOVA dengan menambahkan variabel Age sebagai kovariat agar pengaruh usia dapat dikontrol dalam analisis.
summary(mancova_model, test = "Pillai")## Df Pillai approx F num Df den Df Pr(>F)
## Category 4 1.38737 30.5869 40 2304 < 2.2e-16 ***
## Sex 1 0.16018 10.9291 10 573 < 2.2e-16 ***
## Age 1 0.09328 5.8951 10 573 1.655e-08 ***
## Residuals 582
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1summary(mancova_model, test = "Wilks")## Df Wilks approx F num Df den Df Pr(>F)
## Category 4 0.12200 40.317 40 2174.6 < 2.2e-16 ***
## Sex 1 0.83982 10.929 10 573.0 < 2.2e-16 ***
## Age 1 0.90672 5.895 10 573.0 1.655e-08 ***
## Residuals 582
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1summary.aov(mancova_model)## Response ALB :
## Df Sum Sq Mean Sq F value Pr(>F)
## Category 4 4587.0 1146.75 47.1312 < 2.2e-16 ***
## Sex 1 585.1 585.08 24.0466 1.222e-06 ***
## Age 1 187.8 187.80 7.7187 0.005641 **
## Residuals 582 14160.7 24.33
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Response ALP :
## Df Sum Sq Mean Sq F value Pr(>F)
## Category 4 51143 12785.7 21.9781 < 2.2e-16 ***
## Sex 1 9 9.4 0.0162 0.898902
## Age 1 5349 5349.1 9.1950 0.002535 **
## Residuals 582 338577 581.7
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Response ALT :
## Df Sum Sq Mean Sq F value Pr(>F)
## Category 4 49478 12369.4 36.2725 < 2.2e-16 ***
## Sex 1 6962 6962.2 20.4162 7.548e-06 ***
## Age 1 1030 1029.7 3.0195 0.0828 .
## Residuals 582 198469 341.0
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Response AST :
## Df Sum Sq Mean Sq F value Pr(>F)
## Category 4 276232 69058 113.4684 < 2.2e-16 ***
## Sex 1 4638 4638 7.6201 0.005954 **
## Age 1 96 96 0.1583 0.690846
## Residuals 582 354210 609
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Response BIL :
## Df Sum Sq Mean Sq F value Pr(>F)
## Category 4 62883 15720.8 80.5282 < 2.2e-16 ***
## Sex 1 1468 1468.0 7.5196 0.006291 **
## Age 1 187 187.4 0.9599 0.327609
## Residuals 582 113619 195.2
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Response CHE :
## Df Sum Sq Mean Sq F value Pr(>F)
## Category 4 590.54 147.636 40.3227 < 2.2e-16 ***
## Sex 1 101.35 101.346 27.6798 2.015e-07 ***
## Age 1 0.07 0.073 0.0198 0.8881
## Residuals 582 2130.91 3.661
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Response CHOL :
## Df Sum Sq Mean Sq F value Pr(>F)
## Category 4 69.41 17.3531 15.3752 5.813e-12 ***
## Sex 1 0.14 0.1431 0.1268 0.7219
## Age 1 23.00 23.0001 20.3785 7.693e-06 ***
## Residuals 582 656.87 1.1286
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Response CREA :
## Df Sum Sq Mean Sq F value Pr(>F)
## Category 4 138855 34714 15.1990 7.916e-12 ***
## Sex 1 36917 36917 16.1635 6.575e-05 ***
## Age 1 6235 6235 2.7299 0.09903 .
## Residuals 582 1329262 2284
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Response GGT :
## Df Sum Sq Mean Sq F value Pr(>F)
## Category 4 455756 113939 53.0551 < 2.2e-16 ***
## Sex 1 13428 13428 6.2527 0.012674 *
## Age 1 14801 14801 6.8921 0.008885 **
## Residuals 582 1249880 2148
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Response PROT :
## Df Sum Sq Mean Sq F value Pr(>F)
## Category 4 2705.6 676.41 28.3173 < 2e-16 ***
## Sex 1 59.4 59.35 2.4847 0.11550
## Age 1 156.0 155.97 6.5295 0.01086 *
## Residuals 582 13902.1 23.89
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1Hasilnya kemudian dianalisis kembali menggunakan Pillai Trace untuk melihat pengaruh setelah mempertimbangkan kovariat, selanjutnya digunakan Wilks Lambda untuk memperkuat hasil pengujian multivariat, sebagai angkah akhir, ditinjau kembali pengaruh masing-masing variabel dependen setelah model mencakup kovariat.