cat('ANOVA DUA JALUR')
## ANOVA DUA JALUR
tinggi<-c(25,28,26,27,20,22,21,23,30,32,31,29,
         30,32,31,33,25,27,26,28,35,37,36,34,
         28,30,29,31,23,25,24,26,33,35,34,32)
pupuk <- factor(rep(c("Organik", "Anorganik", "Campuran"), each=12))

tanah <- factor(rep(rep(c("Lempung", "Pasir", "Humus"), each=4),3))

model<-aov(tinggi~pupuk*tanah)
summary(model)
##             Df Sum Sq Mean Sq F value   Pr(>F)    
## pupuk        2    152   76.00    45.6 2.20e-09 ***
## tanah        2    488  244.00   146.4 3.22e-15 ***
## pupuk:tanah  4      0    0.00     0.0        1    
## Residuals   27     45    1.67                     
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
interaction.plot(pupuk,tanah,tinggi,type='b',
                 col=c('red','blue'))

cat('MANOVA\n')
## MANOVA
library(heplots)

pendapatan <- c(3,4,3.5,4.5,3.8,
                6,7,6.5,7.5,6.8,
                10,12,11,13,11.5)
kepuasan <- c(6,7,6.5,7,6.8,
              8,9,8.5,9,8.8,
              9,10,9.5,10,9.8)
pendidikan <- factor(rep(c('SMA','S1','S2'), each=5))

data <- data.frame(
  pendidikan,
  pendapatan,
  kepuasan
)
Y <- cbind(pendapatan, kepuasan)

# Model MANOVA
model_manova <- manova(Y ~ pendidikan)

# Matriks B (Between Group)
rata2 <- aggregate(cbind(pendapatan, kepuasan) ~ pendidikan, data = data, mean)
rata2_total <- colMeans(data[, c("pendapatan", "kepuasan")])

B <- matrix(0, 2, 2)
for (i in 1:3) {
  selisih <- as.matrix(rata2[i, 2:3] - rata2_total)
  B <- B + 5 * (t(selisih) %*% selisih)
}
cat('Matriks B:\n')
## Matriks B:
print(B)
##            pendapatan kepuasan
## pendapatan    152.292 56.60000
## kepuasan       56.600 23.33333
# Matriks W (Within Group)
W <- crossprod(residuals(model_manova))
cat('Matriks W:\n')
## Matriks W:
print(W)
##            pendapatan kepuasan
## pendapatan      7.504    3.514
## kepuasan        3.514    2.136
# Wilks' Lambda secara manual
lambda <- det(W) / det(B + W)
cat("Wilks' Lambda =", lambda, "\n")
## Wilks' Lambda = 0.008067319
# Uji MANOVA Wilks
cat("\nHasil MANOVA - Wilks:\n")
## 
## Hasil MANOVA - Wilks:
summary(model_manova, test = 'Wilks')
##            Df     Wilks approx F num Df den Df   Pr(>F)    
## pendidikan  2 0.0080673   55.735      4     22 3.38e-11 ***
## Residuals  12                                              
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
# Uji MANOVA Pillai
cat("\nHasil MANOVA - Pillai:\n")
## 
## Hasil MANOVA - Pillai:
summary(model_manova, test = 'Pillai')
##            Df Pillai approx F num Df den Df    Pr(>F)    
## pendidikan  2  1.759   43.784      4     24 1.085e-10 ***
## Residuals  12                                            
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
# Uji Box's M
cat("\nUji Box's M:\n")
## 
## Uji Box's M:
boxM(Y, pendidikan)
## 
##  Box's M-test for Homogeneity of Covariance Matrices 
## 
## data:  Y by pendidikan 
## Chi-Sq (approx.) = 7.1899, df = 6, p-value = 0.3036