data_kedelai <- data.frame(
  perlakuan = factor(rep(c("H0", "H1", "H2", "H3", "H4", "H5"), each = 4)),
  produksi = c(
    8.0, 8.1, 7.5, 7.7,
    8.3, 8.2, 8.3, 7.9,
    8.9, 8.1, 8.3, 8.0,
    9.3, 9.0, 8.2, 8.7,
    9.7, 9.0, 8.8, 9.0,
    9.5, 8.9, 8.5, 8.9
  )
)

head(data_kedelai)
##   perlakuan produksi
## 1        H0      8.0
## 2        H0      8.1
## 3        H0      7.5
## 4        H0      7.7
## 5        H1      8.3
## 6        H1      8.2
aggregate(produksi ~ perlakuan,
          data = data_kedelai,
          FUN = mean)
##   perlakuan produksi
## 1        H0    7.825
## 2        H1    8.175
## 3        H2    8.325
## 4        H3    8.800
## 5        H4    9.125
## 6        H5    8.950
model <- aov(produksi ~ perlakuan, data = data_kedelai)

summary(model)
##             Df Sum Sq Mean Sq F value   Pr(>F)    
## perlakuan    5  5.073  1.0147   7.424 0.000618 ***
## Residuals   18  2.460  0.1367                     
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
TukeyHSD(model)
##   Tukey multiple comparisons of means
##     95% family-wise confidence level
## 
## Fit: aov(formula = produksi ~ perlakuan, data = data_kedelai)
## 
## $perlakuan
##         diff         lwr       upr     p adj
## H1-H0  0.350 -0.48075883 1.1807588 0.7605222
## H2-H0  0.500 -0.33075883 1.3307588 0.4262251
## H3-H0  0.975  0.14424117 1.8057588 0.0162551
## H4-H0  1.300  0.46924117 2.1307588 0.0011778
## H5-H0  1.125  0.29424117 1.9557588 0.0048547
## H2-H1  0.150 -0.68075883 0.9807588 0.9915758
## H3-H1  0.625 -0.20575883 1.4557588 0.2109303
## H4-H1  0.950  0.11924117 1.7807588 0.0198257
## H5-H1  0.775 -0.05575883 1.6057588 0.0756965
## H3-H2  0.475 -0.35575883 1.3057588 0.4799701
## H4-H2  0.800 -0.03075883 1.6307588 0.0629634
## H5-H2  0.625 -0.20575883 1.4557588 0.2109303
## H4-H3  0.325 -0.50575883 1.1557588 0.8102934
## H5-H3  0.150 -0.68075883 0.9807588 0.9915758
## H5-H4 -0.175 -1.00575883 0.6557588 0.9831831
boxplot(
  produksi ~ perlakuan,
  data = data_kedelai,
  xlab = "Konsentrasi Hormon",
  ylab = "Produksi Kedelai",
  main = "Pengaruh Konsentrasi Hormon terhadap Produksi Kedelai"
)

shapiro.test(residuals(model))
## 
##  Shapiro-Wilk normality test
## 
## data:  residuals(model)
## W = 0.95294, p-value = 0.3133
bartlett.test(produksi ~ perlakuan,
              data = data_kedelai)
## 
##  Bartlett test of homogeneity of variances
## 
## data:  produksi by perlakuan
## Bartlett's K-squared = 2.4671, df = 5, p-value = 0.7814
par(mfrow = c(2, 2))
plot(model)