One-Way Manova

Berikut merupakan data ukuran bunga Iris yang dibagi menjadi 3 species, yaitu Iris Setosa, Iris Virginica, dan Iris Versicolor. Data ini diperoleh dari website kaggle dengan tautannya sebagai berikut : https://www.kaggle.com/datasets/uciml/iris

library(openxlsx)
data1 <- read.xlsx(file.choose())
head(data1)
##   SepalLengthCm SepalWidthCm PetalLengthCm PetalWidthCm         Species
## 1           6.7          3.3           5.7          2.1  Iris-virginica
## 2           4.7          3.2           1.6          0.2     Iris-setosa
## 3           5.5          3.5           1.3          0.2     Iris-setosa
## 4           6.9          3.1           4.9          1.5 Iris-versicolor
## 5           5.8          2.7           5.1          1.9  Iris-virginica
## 6           5.8          2.7           4.1          1.0 Iris-versicolor

Selanjutnya akan dianalisis apakah perbedaan species pada bunga Iris berpengaruh terhadap panjang sepal, lebar sepal, panjang petal, dan lebar petal nya.

Species : Faktor (Setosa, Virginica, dan Versicolor)
SepalLengthCm : Panjang sepal dengan satuan cm
SepalWidthCm : Lebar sepal dengan satuan cm
PetalLengthCm : Panjang petal dengan satuan cm
PetalWidthCm : Lebar petal dengan satuan cm

Uji Normalitas

Akan dilakukan pengujian asumsi normalitas multivariat menggunakan Mardia’s test.

Hipotesis
H0 : Data berdistribusi normal
H1 : Data tidak berdistribusi normal

Taraf signifikansi
⍺ = 5% = 0.05

Statistik Uji

x1 <- data1[,1]
x2 <- data1[,2]
x3 <- data1[,3]
x4 <- data1[,4]
data1_fix <- data.frame(x1=x1, x2=x2, x3=x3, x4=x4)
library(MVN)
test = mvn(data1_fix, mvnTest = "mardia", univariateTest = "SW", multivariatePlot = "qq")

test
## $multivariateNormality
##              Test          Statistic              p value Result
## 1 Mardia Skewness   66.4979524394953 6.72212639729624e-07     NO
## 2 Mardia Kurtosis -0.196550178829772    0.844179561490338    YES
## 3             MVN               <NA>                 <NA>     NO
## 
## $univariateNormality
##           Test  Variable Statistic   p value Normality
## 1 Shapiro-Wilk    x1        0.9761  0.0102      NO    
## 2 Shapiro-Wilk    x2        0.9838  0.0752      YES   
## 3 Shapiro-Wilk    x3        0.8764  <0.001      NO    
## 4 Shapiro-Wilk    x4        0.9026  <0.001      NO    
## 
## $Descriptives
##      n     Mean   Std.Dev Median Min Max 25th 75th       Skew   Kurtosis
## x1 150 5.843333 0.8280661   5.80 4.3 7.9  5.1  6.4  0.3086407 -0.6058125
## x2 150 3.054000 0.4335943   3.00 2.0 4.4  2.8  3.3  0.3274013  0.1983681
## x3 150 3.758667 1.7644204   4.35 1.0 6.9  1.6  5.1 -0.2689994 -1.4166832
## x4 150 1.198667 0.7631607   1.30 0.1 2.5  0.3  1.8 -0.1029060 -1.3573684

Dari perhitungan didapatkan P-Value :
Mardia Skewness = 6.772e-07
Mardia Kurtosis = 0.844179

Kriteria Uji
Terima H0 apabila P-Value > ⍺
Tolak H0 dalam hal lainnya

Keputusan
Karena P-Value Mardia Skewness (6.772e-07) < ⍺ (0.05), maka Tolak H0

Kesimpulan
Dengan taraf signifikansi 5%, dapat disimpulkan bahwa data tidak memenuhi distribusi Normal.

Akan dilakukan pemotongan data untuk mengecek kembali apakah data berdistribusi Normal

data1_p <- data1_fix[1:50,]
dim(data1_p)
## [1] 50  4
x1_p <- data1[1:50,1]
x2_p <- data1[1:50,2]
x3_p <- data1[1:50,3]
x4_p <- data1[1:50,4]

library(MVN)
test = mvn(data1_p, mvnTest = "mardia", univariateTest = "SW", multivariatePlot = "qq")

test
## $multivariateNormality
##              Test         Statistic            p value Result
## 1 Mardia Skewness  28.8225304887005 0.0913130336328653    YES
## 2 Mardia Kurtosis 0.212893083925202  0.831410356961591    YES
## 3             MVN              <NA>               <NA>    YES
## 
## $univariateNormality
##           Test  Variable Statistic   p value Normality
## 1 Shapiro-Wilk    x1        0.9628    0.1160    YES   
## 2 Shapiro-Wilk    x2        0.9766    0.4208    YES   
## 3 Shapiro-Wilk    x3        0.8934    0.0003    NO    
## 4 Shapiro-Wilk    x4        0.8970    0.0004    NO    
## 
## $Descriptives
##     n  Mean   Std.Dev Median Min Max  25th  75th       Skew   Kurtosis
## x1 50 5.892 0.8891271   5.75 4.3 7.9 5.200 6.475  0.4918990 -0.5790325
## x2 50 3.062 0.4040055   3.00 2.2 4.0 2.825 3.375  0.1897102 -0.4995008
## x3 50 3.814 1.8050960   4.25 1.0 6.9 1.600 5.075 -0.2498570 -1.3393903
## x4 50 1.172 0.7502489   1.35 0.1 2.5 0.200 1.750 -0.1520318 -1.3495766

Didapat nilai P-Value :
Mardia Skewness = 0.09131
Mardia Kurtosis = 0.83141

Keputusan
Karena P-Value Skewness dan Kurtosis nya sudah > ⍺, maka H0 diterima

Kesimpulan
Dengan taraf signifikansi 5% dapat disimpulkan bahwa data berdistribusi Normal Multivariat

Uji Homogenitas

Hipotesis
H0 : s1=s2=s3, matriks kovarians grup adalah sama.
H1 : Minimal ada satu matriks kovarians grup (sk) yang berbeda

Taraf Signifikansi
⍺ : 5%

Statistik Uji

library(MASS)
library(biotools)
## ---
## biotools version 4.2
grup <- data1$Species
head(grup)
## [1] "Iris-virginica"  "Iris-setosa"     "Iris-setosa"     "Iris-versicolor"
## [5] "Iris-virginica"  "Iris-versicolor"
boxM(data = data1_fix, grouping = grup)
## 
##  Box's M-test for Homogeneity of Covariance Matrices
## 
## data:  data1_fix
## Chi-Sq (approx.) = 139.24, df = 20, p-value < 2.2e-16

Dari hasil tersebut didaptkan :
P-Value = 2.2e-16

Kriteria Uji
Terima H0 : P-Value > alpha, tolak dalam hal lainnya

Keputusan
Karena P-Value < alpha, maka tolak H0

Kesimpulan
Dengan taraf signifikansi 5% didapatkan hasil bahwa H0 ditolak yang menandakan bahwa minimal terdapat satu matriks kovarians grup (sk) yang berbeda

Untuk pembelajaran data akan diasumsikan memiliki matriks kovarians grup yang sama

MANOVA

Hipotesis
H0 : μ1 = μ2 = 0 (Faktor Species tidak berpengaruh)
H1 : Terdapat minimal satu μi tidak sama dengan 0, i = 1,2 (Faktor Species berpengaruh)

Taraf Signifikansi
⍺ = 5%

Statistik Uji

dataA <- data1[1:50,]
owm = manova(cbind(dataA$SepalLengthCm, dataA$SepalWidthCm, dataA$PetalLengthCm, dataA$PetalWidthCm)~dataA$Species)
summary(owm)
##               Df Pillai approx F num Df den Df   Pr(>F)    
## dataA$Species  2 1.2438   18.506      8     90 4.24e-16 ***
## Residuals     47                                           
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Dari hasil tersebut didapatkan :
P-Value = 4.24e-16

Kriteria Uji
Terima H0 : P-Value > alpha, tolak dalam hal lainnya

Keputusan
Karena P-Value < alpha, maka tolak H0

Kesimpulan
Dengan taraf signifikansi 5% didapatkan hasil bahwa H0 ditolak yang menandakan bahwa terdapat minimal satu μi tidak sama dengan 0, i = 1,2 (Faktor species berpengaruh)

Uji Lanjut

summary.aov(owm)
##  Response 1 :
##               Df Sum Sq Mean Sq F value    Pr(>F)    
## dataA$Species  2 24.890 12.4452  42.243 3.166e-11 ***
## Residuals     47 13.847  0.2946                      
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
##  Response 2 :
##               Df Sum Sq Mean Sq F value   Pr(>F)    
## dataA$Species  2 3.7752 1.88758   21.01 3.03e-07 ***
## Residuals     47 4.2226 0.08984                     
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
##  Response 3 :
##               Df  Sum Sq Mean Sq F value    Pr(>F)    
## dataA$Species  2 151.363  75.681  428.68 < 2.2e-16 ***
## Residuals     47   8.298   0.177                      
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
##  Response 4 :
##               Df  Sum Sq Mean Sq F value    Pr(>F)    
## dataA$Species  2 26.2160  13.108  451.39 < 2.2e-16 ***
## Residuals     47  1.3648   0.029                      
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Dari hasil tersebut didapatkan :
Semua P-Valaue < alpha

Kesimpulan
Species berpengaruh secara signifikan terhadap Panjang Sepal, Lebar Sepal, Panjang Petal, dan Lebar Petal.

Two-Way Manova

Data berikut merupakan data kesehatan individu dengan variabel independennya Jenis Kelamin serta Perokok/bukan
tautan data : https://www.kaggle.com/datasets/khan1803115/hypertension-risk-model-main

data2 <- read.xlsx(file.choose())
head(data2)
##   male age currentSmoker cigsPerDay BPMeds diabetes totChol sysBP diaBP   BMI
## 1    1  39             0          0      0        0     195 106.0    70 26.97
## 2    0  46             0          0      0        0     250 121.0    81 28.73
## 3    1  48             1         20      0        0     245 127.5    80 25.34
## 4    0  61             1         30      0        0     225 150.0    95 28.58
## 5    0  46             1         23      0        0     285 130.0    84 23.10
## 6    0  43             0          0      0        0     228 180.0   110 30.30
##   heartRate glucose Risk
## 1        80      77    0
## 2        95      76    0
## 3        75      70    0
## 4        65     103    1
## 5        85      85    0
## 6        77      99    1

Selanjutnya akan dianalisis apakah terdapat pengaruh jenis kelamin atau perokok/bukan terhadap variabel-variabel dependennya

male (gender)
age (age of the individual)
currentSmoker (smoking status)
cigsPerDay (number of cigarettes smoked per day)
BPMeds (blood pressure medication usage)
diabetes (diabetes status)
totChol (total cholesterol level)
sysBP (systolic blood pressure)
diaBP (diastolic blood pressure)
BMI (body mass index)
heartRate (heart rate)
glucose (glucose level)
Risk (hypertension risk status)

Uji Normalitas

Akan dilakukan pengujian asumsi normalitas multivariat menggunakan Mardia’s test

Hipotesis
H0 : Data berdistribusi normal
H1 : Data tidak berdistribusi normal

Taraf Signifikansi
⍺ = 5% = 0.05

Statistik UJi

x1 <- data2[,7]
x2 <- data2[,8]
x3 <- data2[,9]
x4 <- data2[,10]
x5 <- data2[,11]
x6 <- data2[,12]
data2_fix <- data.frame(x1=x1, x2=x2, x3=x3, x4=x4, x5=x5, x6=x6)
library(MVN)
test = mvn(data2_fix, mvnTest = "mardia", univariateTest = "SW", multivariatePlot = "qq")

test
## $multivariateNormality
##              Test        Statistic p value Result
## 1 Mardia Skewness 8640.56907024134       0     NO
## 2 Mardia Kurtosis 127.263796387722       0     NO
## 3             MVN             <NA>    <NA>     NO
## 
## $univariateNormality
##           Test  Variable Statistic   p value Normality
## 1 Shapiro-Wilk    x1        0.9710  <0.001      NO    
## 2 Shapiro-Wilk    x2        0.9244  <0.001      NO    
## 3 Shapiro-Wilk    x3        0.9729  <0.001      NO    
## 4 Shapiro-Wilk    x4        0.9658  <0.001      NO    
## 5 Shapiro-Wilk    x5        0.9780  <0.001      NO    
## 6 Shapiro-Wilk    x6        0.6822  <0.001      NO    
## 
## $Descriptives
##       n      Mean   Std.Dev  Median    Min   Max   25th     75th      Skew
## x1 1842 236.94408 44.005896 235.000 135.00 600.0 207.00 263.0000 0.7892911
## x2 1842 132.06026 21.953612 128.000  85.50 295.0 117.00 142.8750 1.3020205
## x3 1842  82.72883 11.916952  82.000  51.00 142.5  74.00  90.0000 0.7309116
## x4 1842  25.71477  4.000058  25.315  15.54  45.8  23.03  27.9875 0.8434938
## x5 1842  75.17318 11.757918  75.000  44.00 140.0  67.00  81.0000 0.5793917
## x6 1842  81.45385 19.635524  78.000  40.00 325.0  72.00  87.0000 4.7328360
##      Kurtosis
## x1  3.0232043
## x2  3.2324724
## x3  1.3441141
## x4  1.7823184
## x5  0.8825237
## x6 39.9435201

Dari perhitungan didapatkan P-Value :
Mardia Skewness = 0
Mardia Kurtosis = 0

Kriteria Uji
Terima H0 apabila P-Value > ⍺
Tolak H0 dalam hal lainnya

Keputusan
Karena P-Value Skewness dan Kurtosis nya < ⍺ (0.05), maka Tolak H0

Kesimpulan
Dengan taraf signifikansi 5%, dapat disimpulkan bahwa data tidak memenuhi distribusi Normal.

Akan dilakukan pemotongan data untuk mengecek kembali apakah data berdistribusi Normal

data2_p <- data2_fix[1:30,]
dim(data2_p)
## [1] 30  6
x1_p <- data2[1:30,7]
x2_p <- data2[1:30,8]
x3_p <- data2[1:30,9]
x4_p <- data2[1:30,10]
x5_p <- data2[1:30,11]
x6_p <- data2[1:30,12]

library(MVN)
test = mvn(data2_p, mvnTest = "mardia", univariateTest = "SW", multivariatePlot = "qq")

test
## $multivariateNormality
##              Test          Statistic           p value Result
## 1 Mardia Skewness   52.3204534855503 0.614935944281329    YES
## 2 Mardia Kurtosis -0.639881846692902 0.522249416652925    YES
## 3             MVN               <NA>              <NA>    YES
## 
## $univariateNormality
##           Test  Variable Statistic   p value Normality
## 1 Shapiro-Wilk    x1        0.9528    0.2006    YES   
## 2 Shapiro-Wilk    x2        0.9483    0.1526    YES   
## 3 Shapiro-Wilk    x3        0.9195    0.0260    NO    
## 4 Shapiro-Wilk    x4        0.9496    0.1651    YES   
## 5 Shapiro-Wilk    x5        0.9523    0.1942    YES   
## 6 Shapiro-Wilk    x6        0.9406    0.0943    YES   
## 
## $Descriptives
##     n      Mean  Std.Dev Median    Min    Max    25th     75th      Skew
## x1 30 245.63333 38.02402 239.50 190.00 332.00 222.000 271.5000 0.4379441
## x2 30 132.51667 20.41445 132.00 100.00 182.00 121.250 141.1250 0.5050106
## x3 30  85.38333 12.29009  85.25  68.00 121.00  78.000  90.0000 0.8441608
## x4 30  26.19733  3.77460  26.20  20.77  34.17  23.135  28.4725 0.4379334
## x5 30  76.63333 10.47323  75.00  60.00  98.00  71.250  83.7500 0.3604811
## x6 30  78.96667 12.05872  77.50  61.00 113.00  70.500  85.0000 0.7946572
##       Kurtosis
## x1 -0.82249779
## x2  0.09310448
## x3  0.79851558
## x4 -0.89439032
## x5 -0.72747921
## x6  0.47630107

Didapat nilai P-Value :
Mardia Skewness = 0.61493
Mardia Kurtosis = 0.52224

Keputusan
Karena P-Value Skewness dan Kurtosis nya sudah > ⍺, maka H0 diterima

Kesimpulan
Dengan taraf signifikansi 5% dapat disimpulkan bahwa data berdistribusi Normal Multivariat

Uji Homogenitas

Hipotesis
H0 : s1=s2=s3, matriks kovarians grup adalah sama. H1 : Minimal ada satu matriks kovarians grup (sk) yang berbeda

Taraf Signifikansi
⍺ : 5%

Statistik Uji

# Gender 
head(data2[1:30,1])
## [1] 1 0 1 0 0 0
boxM(data = data2_p, grouping = data2[1:30,1])
## 
##  Box's M-test for Homogeneity of Covariance Matrices
## 
## data:  data2_p
## Chi-Sq (approx.) = 30.368, df = 21, p-value = 0.08484
# Smoke
head(data2[1:30,3])
## [1] 0 0 1 1 1 0
boxM(data = data2_p, grouping = data2[1:30,3])
## 
##  Box's M-test for Homogeneity of Covariance Matrices
## 
## data:  data2_p
## Chi-Sq (approx.) = 23.758, df = 21, p-value = 0.3049

Dari hasil tersebut didaptkan :
P-Value Gender = 0.08484
P-Value Smoke = 0.3049

Kriteria Uji
Terima H0 : P-Value > alpha, tolak dalam hal lainnya

Keputusan
Karena P-Value > alpha, maka terima H0

Kesimpulan
Dengan taraf signifikansi 5% didapatkan hasil bahwa H0 diterima yang menandakan bahwa data tersebut memiliki matriks kovarians yang sama.

Manova

Hipotesis
H0’ : α1 = α2 = 0 (Faktor gender tidak berpengaruh)
H1’ : Setidaknya ada satu αi yang tidak sama dengan 0 (Faktor gender berpengaruh)

H0’’ : β1 = β2 = 0 (Faktor perokok tidak berpengaruh)
H1’’ : Setidaknya ada satu βj yang tidak sama dengan 0 (Faktor perokok berpengaruh)

H0’’’ : αβij = 0, i =1,2 j=1,2 (Interaksi antara gender dan perokok tidak berpengaruh)
H1’’’ : Setidaknya ada satu αβij yang tidak sama dengan 0 (Interaksi antara gender dan perokok berpengaruh)

Taraf Signifikansi
⍺ = 5%

gender <- as.factor(data2$male[1:30])
perokok <- as.factor(data2$currentSmoker[1:30])
manova <- manova(cbind(x1_p, x2_p, x3_p, x4_p, x5_p, x6_p) ~ gender * perokok, data = data2_p)

summary(manova)
##                Df  Pillai approx F num Df den Df  Pr(>F)  
## gender          1 0.18669   0.8034      6     21 0.57839  
## perokok         1 0.50309   3.5436      6     21 0.01395 *
## gender:perokok  1 0.24528   1.1375      6     21 0.37546  
## Residuals      26                                         
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Dari hasil perhitungan didapatkan P-Value :
Gender : 0.57839
Perokok : 0.01395
Gender:Perokok : 0.37546

Kriteria Uji
Terima H0 : P-Value > alpha, tolak dalam hal lainnya

Keputusan
Karena P-Value perokok < alpha, maka tolak H0

Kesimpulan
Dengan taraf signifikansi 5% didapatkan hasil bahwa untuk faktor perokok H0 ditolak, hasil ini menandakan bahwa setidaknya terdapat minimal satu pengaruh dari faktor perokok terhadap variabel dependennya.

Untuk melihat variabel mana yang dipengaruhi, maka akan dilakukan Uji Lanjut

Uji Lanjut

summary.aov(manova)
##  Response x1_p :
##                Df Sum Sq Mean Sq F value Pr(>F)
## gender          1     10    9.80  0.0067 0.9356
## perokok         1    587  586.62  0.3988 0.5332
## gender:perokok  1   3090 3089.85  2.1007 0.1592
## Residuals      26  38243 1470.87               
## 
##  Response x2_p :
##                Df  Sum Sq Mean Sq F value Pr(>F)
## gender          1   290.1  290.07  0.7457 0.3958
## perokok         1  1126.0 1125.95  2.8944 0.1008
## gender:perokok  1   555.5  555.46  1.4279 0.2429
## Residuals      26 10114.3  389.01               
## 
##  Response x3_p :
##                Df Sum Sq Mean Sq F value Pr(>F)
## gender          1    2.3   2.335  0.0146 0.9046
## perokok         1   70.8  70.755  0.4435 0.5113
## gender:perokok  1  159.6 159.588  1.0004 0.3264
## Residuals      26 4147.7 159.526               
## 
##  Response x4_p :
##                Df  Sum Sq Mean Sq F value   Pr(>F)   
## gender          1  10.220  10.220  1.0546 0.313906   
## perokok         1 126.700 126.700 13.0747 0.001262 **
## gender:perokok  1  24.309  24.309  2.5085 0.125321   
## Residuals      26 251.952   9.690                    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
##  Response x5_p :
##                Df  Sum Sq Mean Sq F value Pr(>F)
## gender          1    1.80   1.800  0.0151 0.9032
## perokok         1    3.25   3.253  0.0273 0.8701
## gender:perokok  1   72.85  72.847  0.6104 0.4417
## Residuals      26 3103.07 119.349               
## 
##  Response x6_p :
##                Df Sum Sq Mean Sq F value  Pr(>F)  
## gender          1  544.3  544.27  3.8559 0.06035 .
## perokok         1    2.0    2.04  0.0144 0.90526  
## gender:perokok  1    0.7    0.66  0.0047 0.94607  
## Residuals      26 3670.0  141.15                  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Dari hasil tersebut didapatkan :
P-Value untuk x4 (BMI) < alpha

Ini menandakan bahwa faktor perokok secara signifikan berpengaruh terhadap BMI