One Way ANOVA

Input Data

library(readr)
data1 <- read_csv("skulls.csv")
## Rows: 150 Columns: 6
## ── Column specification ────────────────────────────────────────────────────────
## Delimiter: ","
## chr (1): epoch
## dbl (5): rownames, mb, bh, bl, nh
## 
## ℹ Use `spec()` to retrieve the full column specification for this data.
## ℹ Specify the column types or set `show_col_types = FALSE` to quiet this message.
head(data1)
## # A tibble: 6 × 6
##   rownames epoch      mb    bh    bl    nh
##      <dbl> <chr>   <dbl> <dbl> <dbl> <dbl>
## 1        1 c4000BC   131   138    89    49
## 2        2 c4000BC   125   131    92    48
## 3        3 c4000BC   131   132    99    50
## 4        4 c4000BC   119   132    96    44
## 5        5 c4000BC   136   143   100    54
## 6        6 c4000BC   138   137    89    56

Pada tes one way ANOVA digunakan data mengenai pengukuran tengkorak mesir dari lima zaman yang berbeda, yaitu tahun 4000, 3300, 1850, 200 sebelum masehi dan 150 tahun masehi. Terdapat 150 observasi dan 5 variabel, yaitu epoch (tengkorang mesir dari berbagai zaman), mb (napas maksimum tengkorak), bh (ketinggian basibergmatic tengkorak), bl (panjang basliveolar tengkorak), dan nh (ketinggian nasal tengkorak). Dalam tes one way ANOVA ini ingin dilihat apakah tengkorang mesir dari lima zaman yang berbeda memiliki ukuran-ukuran yang berbeda.

Uji Normalitas Data

library(MVN)
## Warning: package 'MVN' was built under R version 4.4.1
data1_fix <- data1[,3:6]
head(data1_fix)
## # A tibble: 6 × 4
##      mb    bh    bl    nh
##   <dbl> <dbl> <dbl> <dbl>
## 1   131   138    89    49
## 2   125   131    92    48
## 3   131   132    99    50
## 4   119   132    96    44
## 5   136   143   100    54
## 6   138   137    89    56
test = mvn(data1_fix, mvnTest = "mardia", univariateTest = "SW", multivariatePlot = "qq")

test
## $multivariateNormality
##              Test           Statistic           p value Result
## 1 Mardia Skewness    25.2679738909465 0.191370740880957    YES
## 2 Mardia Kurtosis -0.0433467219770312 0.965425147447826    YES
## 3             MVN                <NA>              <NA>    YES
## 
## $univariateNormality
##           Test  Variable Statistic   p value Normality
## 1 Shapiro-Wilk    mb        0.9913    0.4925    YES   
## 2 Shapiro-Wilk    bh        0.9864    0.1503    YES   
## 3 Shapiro-Wilk    bl        0.9940    0.7887    YES   
## 4 Shapiro-Wilk    nh        0.9840    0.0804    YES   
## 
## $Descriptives
##      n      Mean  Std.Dev Median Min Max 25th 75th        Skew   Kurtosis
## mb 150 133.97333 4.890680    134 119 148  131  137 -0.02808971  0.2162337
## bh 150 132.54667 4.939346    133 120 145  129  136 -0.17427582 -0.1485636
## bl 150  96.46000 5.377844     96  81 114   93  100  0.13808882  0.1432677
## nh 150  50.93333 3.207932     51  44  60   49   53  0.08172942 -0.2068898

Untuk melakukan uji One Way ANOVA, perlu dilakukan pengecekan terlebih dahulu apakah data yang digunakan telah berdistribusi normal multivariat. Maka, dilakukan uji normalitas multivariat dengan mardia test. Dengan Hipotesis,

  • H0 : Data berdistribusi normal multivariat

  • H1 : Data tidak berdistribusi normal multivariat

Tolak H0 apabila P-value < alpha (0,05).

Melalui perhitungan syntax R di atas, didapatkan p-value sebagai berikut:

Mardia P-Value Alpha Keputusan
Skewness 0,19137 0,05 Terima H0
Kurtosis 0,96542 0,05 Terima H0

Data berdistribusi normal multivariat apabila P-Value > alpha (0,05). Dari hasil p-value di atas, dapat disimpulkan bahwa kedua nilai p-value > alpha (0,05), maka Terima H0. Sehingga, dapat disimpulkan bahwa data berdistribusi normal multivariat.

Uji Homogenitas

library(biotools)
## Warning: package 'biotools' was built under R version 4.4.1
## Loading required package: MASS
## ---
## biotools version 4.2
grup <- data1$epoch
head(grup)
## [1] "c4000BC" "c4000BC" "c4000BC" "c4000BC" "c4000BC" "c4000BC"
boxM(data = data1_fix, grouping = grup)
## 
##  Box's M-test for Homogeneity of Covariance Matrices
## 
## data:  data1_fix
## Chi-Sq (approx.) = 45.667, df = 40, p-value = 0.2483

Asumsi selanjutnya yang harus terpenuhi pada data adalah homogenitas. Maka, perlu dilakukan pengecekan lagi terhadap kehomegenitas data. Dengan hipotesis,

  • H0 : Matriks kovarians grup adalah sama

  • H1 : Minimal ada satu matriks kovarians grup yang berbeda

Kriteria Uji:

Tolak H0 apabila P-value < alpha (0,05).

Melalui perhitungan syntax R di atas, didapatkan p-value sebagai berikut:

P-Value BoxM Test 0,2483 0,05 Terima H0

Data memenuhi kehomogenitas apabila P-Value > alpha (0,05). Dari hasil p-value di atas, dapat disimpulkan bahwa kedua nilai p-value > alpha (0,05), maka Terima H0. Sehingga, dapat disimpulkan bahwa data memenuhi kehomgenitasan.

Tes One Way MANOVA

owm = manova(cbind(data1$mb, data1$bh, data1$bl, data1$nh)~data1$epoch)
summary(owm)
##              Df  Pillai approx F num Df den Df    Pr(>F)    
## data1$epoch   4 0.35331    3.512     16    580 4.675e-06 ***
## Residuals   145                                             
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Karena sudah memenuhi asumsi normal dan homogenitas, maka bisa dilakukan tes one way ANOVA. Dengan hipotesis sebagai berikut:

  • H0 : Faktor epoch tidak memberikan hasil pengukuran tengkorak mesir dari lima zaman berbeda

  • H1 : Setidaknya ada satu vektor yang dipengaruhi signifikan oleh faktor epoch

Kriteria Uji:

Tolak H0 apabila P-value < alpha (0,05).

Melalui perhitungan syntax R di atas, didapatkan p-value sebagai berikut:

Sehingga, dapat disimpulkan bahwa sudah cukup bukti untuk mengatakan bahwa faktor epoch berpengruh terhadap pengukuran Egyptians skulls. Maka, dapat dilanjutkan ke uji lanjut untuk mengetahui faktor epoch berpengaruh signifikan terhadap masing-masing vektor secara parsial atau tidak.
P-Value One Way ANOVA 4,675e-06 0,05 Tolak H0

Uji Lanjut (Post Hoc Test)

summary.aov(owm)
##  Response 1 :
##              Df  Sum Sq Mean Sq F value    Pr(>F)    
## data1$epoch   4  502.83 125.707  5.9546 0.0001826 ***
## Residuals   145 3061.07  21.111                      
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
##  Response 2 :
##              Df Sum Sq Mean Sq F value  Pr(>F)  
## data1$epoch   4  229.9  57.477  2.4474 0.04897 *
## Residuals   145 3405.3  23.485                  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
##  Response 3 :
##              Df Sum Sq Mean Sq F value    Pr(>F)    
## data1$epoch   4  803.3 200.823  8.3057 4.636e-06 ***
## Residuals   145 3506.0  24.179                      
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
##  Response 4 :
##              Df Sum Sq Mean Sq F value Pr(>F)
## data1$epoch   4   61.2  15.300   1.507 0.2032
## Residuals   145 1472.1  10.153

Karena hasil tes one way MANOVA menunjukan bahwa faktor epoch memberikan hasil pengukuran tengkorak mesir dari lima zaman yang berbeda. Maka, ingin diketahui lebih lanjut faktor epoch memberikan hasil yang berbeda terhadap variabel vektor yang mana. Dilakukan uji lanjut Post-Hoc Test dengan hipotesis sebagai berikut:

  • H0 : Faktor epoch tidak memberikan hasil pengukuran tengkorak mesir dari lima zaman berbeda pada vektor yang diamati

  • H1 : Faktor epoch memberikan hasil pengukuran tengkorak mesir dari lima zaman yang berbeda pada vektor yang diamati

Kriteria Uji:

Tolak H0 apabila P-value < alpha (0,05).

Melalui perhitungan syntax R di atas, didapatkan p-value sebagai berikut:

Dari hasil di atas, didapatkan bahwa epoch (tengkorang mesir dari lima zaman berbeda) memberikan hasil pengukuran yang berbeda terhadap response nafas maksimum, ketinggian basibregmatic, dan panjang basialiveolar pada tengkorak mesir dari lima zaman berbeda. Sedangkan, epoch tidak memberikan hasil pengukuran yang berbeda terhadap response ketinnggian nasal tengkorak mesir dari lima zaman yang berbeda.
Response P-Value alpha Keputusan
1. Nafas Maksimum 0,0001826 0,05 Tolak H0
2. Ketinggian Basibregmatic 0,04897 0,05 Tolak H0
3. Panjang Basialiveolar 4,636e-06 0,05 Tolak H0
4. Ketinggian Nasal 0,2032 0,05 Terima H0

TWO WAY MANOVA

Input Data

library(readr)
data2 <- read_csv("Salaries.csv")
## Rows: 397 Columns: 7
## ── Column specification ────────────────────────────────────────────────────────
## Delimiter: ","
## chr (3): rank, discipline, sex
## dbl (4): rownames, yrs.since.phd, yrs.service, salary
## 
## ℹ Use `spec()` to retrieve the full column specification for this data.
## ℹ Specify the column types or set `show_col_types = FALSE` to quiet this message.
head(data2)
## # A tibble: 6 × 7
##   rownames rank      discipline yrs.since.phd yrs.service sex   salary
##      <dbl> <chr>     <chr>              <dbl>       <dbl> <chr>  <dbl>
## 1        1 Prof      B                     19          18 Male  139750
## 2        2 Prof      B                     20          16 Male  173200
## 3        3 AsstProf  B                      4           3 Male   79750
## 4        4 Prof      B                     45          39 Male  115000
## 5        5 Prof      B                     40          41 Male  141500
## 6        6 AssocProf B                      6           6 Male   97000

Pada tes one way ANOVA digunakan data mengenai professor, professor madya, dan asisten professor dari sebuah kampus di Amerika Serikat pada tahun 2008-2009. Terdapat 397 observasi dan 7 variabel, yaitu rank (tingkatan profesor), discipline (bagian mengajar), yrs.since.phd (tahun semenjak Phd), yrs.of.service (masa kerja), sex (gender), dan salary (gaji). Dalam tes two way ANOVA ini ingin dilihat apakah tingkatan professor dan bagian mengajar dari sebuah kampus di Amerika Serikat pada tahun 2008-2009 berpengaruh secara signifikan terhadap tahun semenjak Phd, masa kerja, dan gaji.

Uji Normalitas Data

## Uji Normalitas Multivariate 
x1 <- data2[,4]
x2 <- data2[,5]
x3 <- data2[,7]

data2_fix <- data.frame(x1=x1, x2=x2, x3=x3)
library(MVN)
test = mvn(data2_fix, mvnTest = "mardia", univariateTest = "SW", multivariatePlot = "qq")

test
## $multivariateNormality
##              Test        Statistic              p value Result
## 1 Mardia Skewness  292.53660804796 5.87259177177525e-57     NO
## 2 Mardia Kurtosis 5.93488545071384 2.94050717108973e-09     NO
## 3             MVN             <NA>                 <NA>     NO
## 
## $univariateNormality
##           Test      Variable Statistic   p value Normality
## 1 Shapiro-Wilk yrs.since.phd    0.9696  <0.001      NO    
## 2 Shapiro-Wilk  yrs.service     0.9418  <0.001      NO    
## 3 Shapiro-Wilk    salary        0.9599  <0.001      NO    
## 
## $Descriptives
##                 n         Mean     Std.Dev Median   Min    Max  25th   75th
## yrs.since.phd 397     22.31486    12.88700     21     1     56    12     32
## yrs.service   397     17.61461    13.00602     16     0     60     7     27
## salary        397 113706.45844 30289.03869 107300 57800 231545 91000 134185
##                    Skew   Kurtosis
## yrs.since.phd 0.2986097 -0.8105942
## yrs.service   0.6456616 -0.3360614
## salary        0.7091778  0.1814838

Untuk melakukan uji Two Way MANOVA, perlu dilakukan pengecekan terlebih dahulu apakah data yang digunakan telah berdistribusi normal multivariat. Maka, dilakukan uji normalitas multivariat dengan mardia test. Dengan Hipotesis,

  • H0 : Data berdistribusi normal multivariat

  • H1 : Data tidak berdistribusi normal multivariat

Tolak H0 apabila P-value < alpha (0,05).

Melalui perhitungan syntax R di atas, didapatkan p-value sebagai berikut:

Mardia P-Value Alpha Keputusan
Skewness 5,87259 0,05 Tolak H0
Kurtosis 2,94051 0,05 Tolak H0

Data berdistribusi normal multivariat apabila Terima H0,P-Value > alpha (0,05).

Dari hasil p-value di atas, dapat disimpulkan bahwa kedua nilai p-value < alpha (0,05), maka Tolak H0. Sehingga, dapat disimpulkan bahwa data tidak berdistribusi normal multivariat. Dengan tujuan pembelajaran, akan dilakukan pemotongan terhadap data yang awalnya 397 menjadi 30 (85:114) agar data menjadi berdistribusi normal multivariat.

Pemotongan Data & Uji Normalitas Data 2

## Potong Data 
data2_p <- data2_fix[85:114,]
dim(data2_p)
## [1] 30  3
x1_p <- data2[85:114,4]
x2_p <- data2[85:114,5]
x3_p <- data2[85:114,7]

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

test
## $multivariateNormality
##              Test         Statistic           p value Result
## 1 Mardia Skewness  18.1303365648699 0.052804447146691    YES
## 2 Mardia Kurtosis 0.453356519733072 0.650292040267642    YES
## 3             MVN              <NA>              <NA>    YES
## 
## $univariateNormality
##           Test      Variable Statistic   p value Normality
## 1 Shapiro-Wilk yrs.since.phd    0.9302    0.0496    NO    
## 2 Shapiro-Wilk  yrs.service     0.9258    0.0381    NO    
## 3 Shapiro-Wilk    salary        0.9461    0.1332    YES   
## 
## $Descriptives
##                n         Mean     Std.Dev   Median   Min    Max     25th
## yrs.since.phd 30     20.30000    11.95437     17.5     2     42    10.25
## yrs.service   30     15.63333    11.08115     14.0     0     38     7.00
## salary        30 112631.73333 26428.28317 106848.5 72500 172272 90202.75
##                   75th      Skew   Kurtosis
## yrs.since.phd     29.5 0.4133458 -1.0711037
## yrs.service       24.5 0.5740170 -0.8426250
## salary        129135.0 0.5074847 -0.6278946

Setelah dilakukan pemotongan data, dilakukan kembali uji normlaitas pada data. Dengan Hipotesis,

  • H0 : Data berdistribusi normal multivariat

  • H1 : Data tidak berdistribusi normal multivariat

Tolak H0 apabila P-value < alpha (0,05).

Melalui perhitungan syntax R di atas, didapatkan p-value sebagai berikut:

Mardia P-Value Alpha Keputusan
Skewness 0,0528 0,05 Terima H0
Kurtosis 0,65029 0,05 Terima H0

Data berdistribusi normal multivariat apabila Terima H0,P-Value > alpha (0,05).

Dari hasil p-value di atas, dapat disimpulkan bahwa kedua nilai p-value > alpha (0,05), maka Tolak H0. Sehingga, dapat disimpulkan bahwa data tidak berdistribusi normal multivariat.

Uji Homogenitas Data

library(biotools)
grup1 <- ifelse(data2[85:114,2] == "Prof",1, ifelse(data2[85:114,2] == "AsstProf",2,3))
head(data2[85:114,2])
## # A tibble: 6 × 1
##   rank     
##   <chr>    
## 1 Prof     
## 2 Prof     
## 3 Prof     
## 4 AsstProf 
## 5 Prof     
## 6 AssocProf
boxM(data = data2_p, grouping = grup1)
## 
##  Box's M-test for Homogeneity of Covariance Matrices
## 
## data:  data2_p
## Chi-Sq (approx.) = 26.216, df = 12, p-value = 0.01
grup2 <- ifelse(data2[85:114,3] == "Prof",1,2)
head(data2[85:114,3])
## # A tibble: 6 × 1
##   discipline
##   <chr>     
## 1 B         
## 2 B         
## 3 B         
## 4 B         
## 5 B         
## 6 B
boxM(data = data2_p, grouping = grup2)
## 
##  Box's M-test for Homogeneity of Covariance Matrices
## 
## data:  data2_p
## Chi-Sq (approx.) = NaN, df = 0, p-value = NA

Asumsi selanjutnya yang harus terpenuhi pada data adalah homogenitas. Maka, perlu dilakukan pengecekan lagi terhadap kehomegenitas data. Dengan hipotesis,

  • H0 : Matriks kovarians grup adalah sama

  • H1 : Minimal ada satu matriks kovarians grup yang berbeda

Kriteria Uji:

Tolak H0 apabila P-value < alpha (0,05).

Melalui perhitungan syntax R di atas, didapatkan p-value sebagai berikut:

P-Value BoxM Test (Rank) 0,01 0,05 Tolak H0
P-Value BoxM Test (Discipline) 0,168 0,05 Terima H0

Data memenuhi kehomogenitas apabila P-Value > alpha (0,05). Dari hasil p-value di atas, dapat disimpulkan bahwa grup rank nilai p-value < alpha (0,05), maka Tolak H0. Sehingga, dapat disimpulkan bahwa data grup rank tidak memenuhi kehomgenitasan. Dari hasil p-value di atas, dapat disimpulkan bahwa grup discipline nilai p-value > alpha (0,05), maka Terima H0. Sehingga, dapat disimpulkan bahwa data grup discipline memenuhi kehomgenitasan. Namun, dengan tujuan pembelajaran untuk grup rank diasumsikan telah memenuhi asumsi.

Tes Two Way ANOVA

## Two Way MANOVA 
rank <- as.factor(data2$rank[85:114])
discipline <- as.factor(data2$discipline[85:114])
x1_p<-unlist(x1_p)
x2_p<-unlist(x2_p)
x3_p<-unlist(x3_p)
manova <- manova(cbind(x1_p, x2_p, x3_p) ~ discipline * rank, data = data2_p)
# Menampilkan hasil
summary(manova)
##                 Df  Pillai approx F num Df den Df    Pr(>F)    
## discipline       1 0.53385   8.3985      3     22 0.0006607 ***
## rank             2 0.75378   4.6372      6     46 0.0009205 ***
## discipline:rank  2 0.06562   0.2601      6     46 0.9525695    
## Residuals       24                                             
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Karena sudah memenuhi asumsi normal dan homogenitas, maka bisa dilakukan tes Two Way MANOVA. Dengan hipotesis sebagai berikut:

Rank:

  • H0 : Faktor rank tidak berpengaruh signifikan terhadap vektor

    H1 : Setidaknya ada satu vektor yang dipengaruhi signifikan oleh faktor rank

Discipline:

  • H0 : Faktor discipline tidak berpengaruh signifikan terhadap vektor

  • H1 : Setidaknya ada satu vektor yang dipengaruhi signifikan oleh faktor discipline

Rank & Discipline:

  • H0 : Interaksi antara faktor rank & discipline tidak berpengaruh signifikan terhadap vektor

  • H1 : Setidaknya ada satu vektor yang dipengaruhi signifikan oleh interaksi antara faktor rank & discipline

Kriteria Uji:

Tolak H0 apabila P-value < alpha (0,05).

Melalui perhitungan syntax R di atas, didapatkan p-value sebagai berikut:

Sehingga, dapat disimpulkan bahwa sudah cukup bukti untuk mengatakan bahwa faktor rank dan discipline berpengruh terhadap vektor satuan. Namun, untuk interaksi antara faktor rank dan discipline tidak berpengaruh terhadap vektor satuan. Maka, dapat dilanjutkan ke uji lanjut untuk mengetahui faktor epoch berpengaruh signifikan terhadap masing-masing vektor secara parsial atau tidak.
P-Value Rank 0,00067 0,05 Tolak H0
P-Value Discipline 0,0009205 0,05 Tolak H0
P-Value Rank & Discipline 0,95257 0,05 Terima H0

Uji Lanjut Post Hoc Test

summary.aov(manova)
##  Response x1_p :
##                 Df  Sum Sq Mean Sq F value    Pr(>F)    
## discipline       1    0.15    0.15  0.0022    0.9634    
## rank             2 2432.88 1216.44 17.4615 2.085e-05 ***
## discipline:rank  2   39.32   19.66  0.2822    0.7566    
## Residuals       24 1671.94   69.66                      
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
##  Response x2_p :
##                 Df  Sum Sq Mean Sq F value    Pr(>F)    
## discipline       1    6.67    6.67  0.1049    0.7488    
## rank             2 1959.32  979.66 15.4212 4.933e-05 ***
## discipline:rank  2   70.33   35.17  0.5535    0.5821    
## Residuals       24 1524.65   63.53                      
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
##  Response x3_p :
##                 Df     Sum Sq    Mean Sq F value    Pr(>F)    
## discipline       1 4871472763 4871472763 19.9437 0.0001616 ***
## rank             2 9513287357 4756643678 19.4736 9.437e-06 ***
## discipline:rank  2    8133878    4066939  0.0166 0.9834993    
## Residuals       24 5862276395  244261516                      
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Kriteria Uji:

Tolak H0 apabila P-value < alpha (0,05).

Response Group P-Value alpha Keputusan
Rank Years since Phd. 2,085e-05 0,05 Tolak H0
Rank Years of service 4,933e-05 0,05 Tolak H0
Rank Salary 9,437e-06 0,05 Tolak H0
Discipline Years since Phd. 0,9643 0,05 Terima H0
Discipline Years of service 0,7488 0,05 Terima H0
Discipline Salary 0,0001616 0,05 Tolak H0
Rank & Discipline Years since Phd. 0,7566 0,05 Terima H0
Rank & Discipline Years of service 0,5821 0,05 Terima H0
Rank & Discipline Salary 0,9835 0,05 Terima H0

Kesimpulan:

  • Faktor rank berpengaruh signifikan terhadap tahun semenjak Phd

  • Faktor rank berpengaruh signifikan terhadap masa kerja

  • Faktor rank berpengaruh signifikan terhadap gaji

  • Faktor discipline tidak berpengaruh signifikan terhadap tahun semenjak Phd

  • Faktor discipline tidak berpengaruh signifikan terhadap masa kerja

  • Faktor discipline berpengaruh signifikan terhadap gaji

  • Faktor interaksi antara rank dan discipline tidak berpengaruh signifikan terhadap tahun semenjak Phd

  • Faktor interaksi rank dan discipline tidak berpengaruh signifikan terhadap masa kerja

  • Faktor interaksi antara rank dan discipline tidak berpengaruh signifikan terhadap gaji