ONE WAY MANOVA ### One Way MANOVA

data <- read.csv(file.choose())
head(data)
##   rownames Active Rest Smoke Sex Exercise Hgt Wgt
## 1        1     97   78     0   1        1  63 119
## 2        2     82   68     1   0        3  70 225
## 3        3     88   62     0   0        3  72 175
## 4        4    106   74     0   0        3  72 170
## 5        5     78   63     0   1        3  67 125
## 6        6    109   65     0   0        3  74 188

Uji Normalita`s

library(MVN)
## Warning: package 'MVN' was built under R version 4.3.3
data1<- data[,7:8]
head(data1)
##   Hgt Wgt
## 1  63 119
## 2  70 225
## 3  72 175
## 4  72 170
## 5  67 125
## 6  74 188
test1 = mvn(data1, mvnTest = "mardia", univariateTest = "SW", multivariatePlot = "qq")

test1
## $multivariateNormality
##              Test        Statistic              p value Result
## 1 Mardia Skewness 42.6121795459257 1.24539244381758e-08     NO
## 2 Mardia Kurtosis 1.57723328287731    0.114741857341793    YES
## 3             MVN             <NA>                 <NA>     NO
## 
## $univariateNormality
##           Test  Variable Statistic   p value Normality
## 1 Shapiro-Wilk    Hgt       0.9820  0.0047      NO    
## 2 Shapiro-Wilk    Wgt       0.9574  <0.001      NO    
## 
## $Descriptives
##       n      Mean   Std.Dev Median Min Max 25th 75th        Skew   Kurtosis
## Hgt 232  68.24569  3.738761     68  60  78   65   71 -0.03945937 -0.6631010
## Wgt 232 157.91810 31.832587    150 102 260  135  175  0.77037561  0.5246982
data1_p<- data[1:40,]
head(data1_p)
##   rownames Active Rest Smoke Sex Exercise Hgt Wgt
## 1        1     97   78     0   1        1  63 119
## 2        2     82   68     1   0        3  70 225
## 3        3     88   62     0   0        3  72 175
## 4        4    106   74     0   0        3  72 170
## 5        5     78   63     0   1        3  67 125
## 6        6    109   65     0   0        3  74 188
test1_p = mvn(data1_p[,7:8], mvnTest = "mardia", univariateTest = "SW", multivariatePlot = "qq")

test1_p
## $multivariateNormality
##              Test          Statistic           p value Result
## 1 Mardia Skewness   4.23864113028207 0.374671605772214    YES
## 2 Mardia Kurtosis 0.0813720009195687 0.935146115282007    YES
## 3             MVN               <NA>              <NA>    YES
## 
## $univariateNormality
##           Test  Variable Statistic   p value Normality
## 1 Shapiro-Wilk    Hgt       0.9553    0.1153    YES   
## 2 Shapiro-Wilk    Wgt       0.9513    0.0841    YES   
## 
## $Descriptives
##      n    Mean   Std.Dev Median Min Max 25th 75th       Skew   Kurtosis
## Hgt 40  68.875  3.567571   69.5  62  75   66   72 -0.0356674 -1.1631583
## Wgt 40 158.050 28.918054  156.5 115 225  135  175  0.3289874 -0.9577016
str(data1_p)
## 'data.frame':    40 obs. of  8 variables:
##  $ rownames: int  1 2 3 4 5 6 7 8 9 10 ...
##  $ Active  : int  97 82 88 106 78 109 66 68 100 70 ...
##  $ Rest    : int  78 68 62 74 63 65 43 65 63 59 ...
##  $ Smoke   : int  0 1 0 0 0 0 0 0 0 0 ...
##  $ Sex     : int  1 0 0 0 1 0 1 0 0 1 ...
##  $ Exercise: int  1 3 3 3 3 3 3 3 1 2 ...
##  $ Hgt     : int  63 70 72 72 67 74 67 70 70 65 ...
##  $ Wgt     : int  119 225 175 170 125 188 140 200 165 115 ...
  ## Uji Homogenitas Multivariate 

library(biotools)
## Warning: package 'biotools' was built under R version 4.3.3
## Loading required package: MASS
## ---
## biotools version 4.2
grup <- (data1_p$Sex)
dim(grup)
## NULL
length(grup)
## [1] 40
boxM(data = data1_p, grouping = grup)
## 
##  Box's M-test for Homogeneity of Covariance Matrices
## 
## data:  data1_p
## Chi-Sq (approx.) = NaN, df = 36, p-value = NA

One Way MANOVA

owm = manova(cbind(data1_p$Wgt, data1_p$Hgt)~data1_p$Sex)
summary(owm)
##             Df  Pillai approx F num Df den Df    Pr(>F)    
## data1_p$Sex  1 0.70349   43.893      2     37 1.709e-10 ***
## Residuals   38                                             
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Uji Lanjut (Post Hoc Test)

summary.aov(owm)
##  Response 1 :
##             Df Sum Sq Mean Sq F value    Pr(>F)    
## data1_p$Sex  1  19017 19017.1  53.148 9.989e-09 ***
## Residuals   38  13597   357.8                      
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
##  Response 2 :
##             Df Sum Sq Mean Sq F value    Pr(>F)    
## data1_p$Sex  1 325.31  325.31  72.263 2.542e-10 ***
## Residuals   38 171.07    4.50                      
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

TWO WAY MANOVA

#Menggunakan data yang sama ## Uji Normalitas Multivariate

library(MVN)
data2<- data[,2:3]
head(data2)
##   Active Rest
## 1     97   78
## 2     82   68
## 3     88   62
## 4    106   74
## 5     78   63
## 6    109   65
test2 = mvn(data2[], mvnTest = "mardia", univariateTest = "SW", multivariatePlot = "qq")

test2
## $multivariateNormality
##              Test        Statistic              p value Result
## 1 Mardia Skewness 56.9297031715916 1.27990072771327e-11     NO
## 2 Mardia Kurtosis 4.44548141755373 8.76952613526605e-06     NO
## 3             MVN             <NA>                 <NA>     NO
## 
## $univariateNormality
##           Test  Variable Statistic   p value Normality
## 1 Shapiro-Wilk  Active      0.9704    0.0001    NO    
## 2 Shapiro-Wilk   Rest       0.9825    0.0057    NO    
## 
## $Descriptives
##          n     Mean   Std.Dev Median Min Max 25th 75th      Skew  Kurtosis
## Active 232 91.29741 18.820234   88.5  51 154   79  102 0.6647994 0.3670902
## Rest   232 68.34914  9.949378   68.0  43 106   62   74 0.5158032 0.7603768

Potong Data

data2_p<- data[1:40,]
head(data2_p)
##   rownames Active Rest Smoke Sex Exercise Hgt Wgt
## 1        1     97   78     0   1        1  63 119
## 2        2     82   68     1   0        3  70 225
## 3        3     88   62     0   0        3  72 175
## 4        4    106   74     0   0        3  72 170
## 5        5     78   63     0   1        3  67 125
## 6        6    109   65     0   0        3  74 188
dim(data2_p)
## [1] 40  8
x1_p <- data2_p[,2]
x2_p <- data2_p[,3]

Uji Normalitas Lagi

library(MVN)
test2_p = mvn(data2_p[,2:3], mvnTest = "mardia", univariateTest = "SW", multivariatePlot = "qq")

test2_p
## $multivariateNormality
##              Test          Statistic           p value Result
## 1 Mardia Skewness   5.55906849547326 0.234585896783726    YES
## 2 Mardia Kurtosis -0.232463096056512 0.816178347237994    YES
## 3             MVN               <NA>              <NA>    YES
## 
## $univariateNormality
##           Test  Variable Statistic   p value Normality
## 1 Shapiro-Wilk  Active      0.9562    0.1241    YES   
## 2 Shapiro-Wilk   Rest       0.9844    0.8445    YES   
## 
## $Descriptives
##         n   Mean  Std.Dev Median Min Max  25th   75th       Skew   Kurtosis
## Active 40 92.975 20.86893   87.5  61 140 77.50 108.25 0.41233985 -0.7494245
## Rest   40 67.700 10.78032   68.0  43  97 61.75  74.00 0.06214363  0.2328359

Two Way MANOVA

Sex <- as.factor(data2_p$Sex)
Sex
##  [1] 1 0 0 0 1 0 1 0 0 1 0 1 1 0 1 1 0 0 0 0 1 1 1 0 0 1 1 0 0 1 1 0 0 0 1 0 1 0
## [39] 1 0
## Levels: 0 1
Smoke <- as.factor(data2_p$Smoke)
Smoke
##  [1] 0 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 0 0 0 0 1 0 0 0 0 0 0 0 1 0 0 0 0 0
## [39] 0 0
## Levels: 0 1
manova <- manova(cbind(x1_p, x2_p) ~ Sex * Smoke, data = data2_p)

Menampilkan hasil

summary(manova)
##           Df   Pillai approx F num Df den Df Pr(>F)
## Sex        1 0.088891  1.75614      2     36 0.1872
## Smoke      1 0.036871  0.68908      2     36 0.5085
## Residuals 37

Uji Lanjut

summary.aov(manova)
##  Response x1_p :
##             Df  Sum Sq Mean Sq F value  Pr(>F)  
## Sex          1  1489.9 1489.91  3.5613 0.06701 .
## Smoke        1    15.8   15.76  0.0377 0.84717  
## Residuals   37 15479.3  418.36                  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
##  Response x2_p :
##             Df Sum Sq Mean Sq F value Pr(>F)
## Sex          1  156.8  156.80  1.3696 0.2494
## Smoke        1  139.5  139.46  1.2181 0.2769
## Residuals   37 4236.1  114.49