Pengertian Uji Lanjut

Uji lanjut (post hoc test) merupakan analisis statistik yang dilakukan setelah uji ANOVA menunjukkan hasil yang signifikan, yaitu ketika hipotesis nol (H₀) ditolak. Penolakan H₀ menunjukkan bahwa terdapat perbedaan rata-rata antar perlakuan, namun uji ANOVA tidak mampu menjelaskan pasangan perlakuan mana yang berbeda secara signifikan.

Oleh karena itu, uji lanjut digunakan untuk mengidentifikasi secara spesifik pasangan perlakuan yang memberikan pengaruh berbeda terhadap rata-rata respon.

Berdasarkan perencanaannya, uji lanjut dibedakan menjadi dua jenis:

1. Uji Terencana (Planned Comparison)

Digunakan apabila perbandingan antar perlakuan telah dirancang sebelum penelitian dilakukan, biasanya menggunakan kontras ortogonal dan Polynomial.

2. Uji Tidak Terencana (Unplanned Comparison)

Digunakan apabila tidak ada perencanaan sebelumnya, sehingga dilakukan eksplorasi terhadap seluruh pasangan perlakuan.
Metode yang digunakan antara lain: - BNT (Beda Nyata Terkecil / LSD) - BNJ (Beda Nyata Jujur / Tukey) - Duncan (DMRT)

Dengan demikian, uji lanjut berfungsi untuk memberikan informasi lebih rinci mengenai letak perbedaan yang tidak dapat dijelaskan oleh ANOVA.

UJI LANJUT TAK TERENCANA

Contoh Aplikasi diRstudio

# Instal dan muat paket
library(agricolae)
## Warning: package 'agricolae' was built under R version 4.3.3
# Data kadar nitrogen dari Contoh Kasus 1
data <- data.frame(
  Perlakuan = rep(c("3DOK1", "3DOK5", "3DOK4", "3DOK7", "3DOK13", "KOMPOSIT"), each = 5),
  Nitrogen = c(
    19.4, 32.6, 27.0, 32.1, 33.0,      # 3DOK1
    17.7, 24.8, 27.9, 25.2, 24.3,      # 3DOK5
    17.0, 19.4, 9.1, 11.9, 15.8,       # 3DOK4
    20.7, 21.0, 20.5, 18.8, 18.6,      # 3DOK7
    14.3, 14.4, 11.8, 11.6, 14.2,      # 3DOK13
    17.3, 19.4, 19.1, 16.9, 20.8       # KOMPOSIT
  )
)

# Cek ringkasan data
head(data)
##   Perlakuan Nitrogen
## 1     3DOK1     19.4
## 2     3DOK1     32.6
## 3     3DOK1     27.0
## 4     3DOK1     32.1
## 5     3DOK1     33.0
## 6     3DOK5     17.7

Melakukan analisis ragam/anova

# Model ANOVA
model <- aov(Nitrogen ~ Perlakuan, data = data)

# Hasil ANOVA
anova(model)
## Analysis of Variance Table
## 
## Response: Nitrogen
##           Df Sum Sq Mean Sq F value    Pr(>F)    
## Perlakuan  5 847.05 169.409   14.37 1.485e-06 ***
## Residuals 24 282.93  11.789                      
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

3. Uji Lanjut

3.1 Uji lanjut BNT

# Uji BNT (Beda Nyata Terkecil)
hasil_BNT <- LSD.test(model, "Perlakuan", group = TRUE)
hasil_BNT
## $statistics
##    MSerror Df     Mean       CV  t.value      LSD
##   11.78867 24 19.88667 17.26515 2.063899 4.481782
## 
## $parameters
##         test p.ajusted    name.t ntr alpha
##   Fisher-LSD      none Perlakuan   6  0.05
## 
## $means
##          Nitrogen      std r       se     LCL     UCL  Min  Max  Q25  Q50  Q75
## 3DOK1       28.82 5.800172 5 1.535491 25.6509 31.9891 19.4 33.0 27.0 32.1 32.6
## 3DOK13      13.26 1.427585 5 1.535491 10.0909 16.4291 11.6 14.4 11.8 14.2 14.3
## 3DOK4       14.64 4.116188 5 1.535491 11.4709 17.8091  9.1 19.4 11.9 15.8 17.0
## 3DOK5       23.98 3.777168 5 1.535491 20.8109 27.1491 17.7 27.9 24.3 24.8 25.2
## 3DOK7       19.92 1.130044 5 1.535491 16.7509 23.0891 18.6 21.0 18.8 20.5 20.7
## KOMPOSIT    18.70 1.601562 5 1.535491 15.5309 21.8691 16.9 20.8 17.3 19.1 19.4
## 
## $comparison
## NULL
## 
## $groups
##          Nitrogen groups
## 3DOK1       28.82      a
## 3DOK5       23.98      b
## 3DOK7       19.92     bc
## KOMPOSIT    18.70     cd
## 3DOK4       14.64     de
## 3DOK13      13.26      e
## 
## attr(,"class")
## [1] "group"

➡️ Huruf yang sama → tidak berbeda nyata ➡️ Huruf berbeda → berbeda nyata 1. 3DOK1 - Tertinggi; berbeda nyata dengan semua perlakuan lain (b–e).

3.2 Uji Lanjut BNJ/Tukey

# Uji BNJ (Tukey HSD)
hasil_BNJ <- HSD.test(model, "Perlakuan", group = TRUE)
hasil_BNJ
## $statistics
##    MSerror Df     Mean       CV      MSD
##   11.78867 24 19.88667 17.26515 6.714167
## 
## $parameters
##    test    name.t ntr StudentizedRange alpha
##   Tukey Perlakuan   6         4.372651  0.05
## 
## $means
##          Nitrogen      std r       se  Min  Max  Q25  Q50  Q75
## 3DOK1       28.82 5.800172 5 1.535491 19.4 33.0 27.0 32.1 32.6
## 3DOK13      13.26 1.427585 5 1.535491 11.6 14.4 11.8 14.2 14.3
## 3DOK4       14.64 4.116188 5 1.535491  9.1 19.4 11.9 15.8 17.0
## 3DOK5       23.98 3.777168 5 1.535491 17.7 27.9 24.3 24.8 25.2
## 3DOK7       19.92 1.130044 5 1.535491 18.6 21.0 18.8 20.5 20.7
## KOMPOSIT    18.70 1.601562 5 1.535491 16.9 20.8 17.3 19.1 19.4
## 
## $comparison
## NULL
## 
## $groups
##          Nitrogen groups
## 3DOK1       28.82      a
## 3DOK5       23.98     ab
## 3DOK7       19.92     bc
## KOMPOSIT    18.70     bc
## 3DOK4       14.64      c
## 3DOK13      13.26      c
## 
## attr(,"class")
## [1] "group"

Artinya:

Perlakuan dengan huruf yang sama → tidak berbeda nyata.

Perlakuan dengan huruf berbeda → berbeda nyata.

  1. 3DOK1- Tertinggi; berbeda nyata dari kelompok c, tapi tidak berbeda nyata dengan 3DOK5 (karena 3DOK5 punya “a” juga).
# Tukey HSD bawaan R
TukeyHSD(model)
##   Tukey multiple comparisons of means
##     95% family-wise confidence level
## 
## Fit: aov(formula = Nitrogen ~ Perlakuan, data = data)
## 
## $Perlakuan
##                   diff          lwr       upr     p adj
## 3DOK13-3DOK1    -15.56 -22.27416704 -8.845833 0.0000029
## 3DOK4-3DOK1     -14.18 -20.89416704 -7.465833 0.0000128
## 3DOK5-3DOK1      -4.84 -11.55416704  1.874167 0.2617111
## 3DOK7-3DOK1      -8.90 -15.61416704 -2.185833 0.0048849
## KOMPOSIT-3DOK1  -10.12 -16.83416704 -3.405833 0.0012341
## 3DOK4-3DOK13      1.38  -5.33416704  8.094167 0.9870716
## 3DOK5-3DOK13     10.72   4.00583296 17.434167 0.0006233
## 3DOK7-3DOK13      6.66  -0.05416704 13.374167 0.0527514
## KOMPOSIT-3DOK13   5.44  -1.27416704 12.154167 0.1621550
## 3DOK5-3DOK4       9.34   2.62583296 16.054167 0.0029837
## 3DOK7-3DOK4       5.28  -1.43416704 11.994167 0.1852490
## KOMPOSIT-3DOK4    4.06  -2.65416704 10.774167 0.4434643
## 3DOK7-3DOK5      -4.06 -10.77416704  2.654167 0.4434643
## KOMPOSIT-3DOK5   -5.28 -11.99416704  1.434167 0.1852490
## KOMPOSIT-3DOK7   -1.22  -7.93416704  5.494167 0.9926132
library(ggplot2)

# Ambil hasil Tukey
tukey_res <- TukeyHSD(model)

# Jadikan data frame
tukey_df <- as.data.frame(tukey_res$Perlakuan)
tukey_df$comparison <- rownames(tukey_df)

# Plot dengan ggplot (lebih teratur dan jelas)
ggplot(tukey_df, aes(x = diff, y = reorder(comparison, diff))) +
  geom_point(size = 3, color = "blue") +
  geom_errorbarh(aes(xmin = lwr, xmax = upr), height = 0.2, color = "gray40") +
  geom_vline(xintercept = 0, linetype = "dashed", color = "red") +
  labs(
    title = "Hasil Uji Tukey HSD (Tersusun Rapi)",
    subtitle = "Interval kepercayaan 95% antar perlakuan",
    x = "Selisih rata-rata antar perlakuan",
    y = "Perbandingan antar perlakuan"
  ) +
  theme_minimal(base_size = 13)
## Warning: `geom_errobarh()` was deprecated in ggplot2 4.0.0.
## ℹ Please use the `orientation` argument of `geom_errorbar()` instead.
## This warning is displayed once every 8 hours.
## Call `lifecycle::last_lifecycle_warnings()` to see where this warning was
## generated.
## `height` was translated to `width`.

3.3 Duncan

Huruf (a, b, c, …) disebut kelompok huruf Duncan (grouping letters). Artinya:

Perlakuan dengan huruf yang sama → tidak berbeda nyata.

Perlakuan dengan huruf berbeda → berbeda nyata.

# Uji Duncan
hasil_Duncan <- duncan.test(model, "Perlakuan", group = TRUE)
hasil_Duncan
## $statistics
##    MSerror Df     Mean       CV
##   11.78867 24 19.88667 17.26515
## 
## $parameters
##     test    name.t ntr alpha
##   Duncan Perlakuan   6  0.05
## 
## $duncan
##      Table CriticalRange
## 2 2.918793      4.481781
## 3 3.065610      4.707218
## 4 3.159874      4.851958
## 5 3.226454      4.954192
## 6 3.276155      5.030507
## 
## $means
##          Nitrogen      std r       se  Min  Max  Q25  Q50  Q75
## 3DOK1       28.82 5.800172 5 1.535491 19.4 33.0 27.0 32.1 32.6
## 3DOK13      13.26 1.427585 5 1.535491 11.6 14.4 11.8 14.2 14.3
## 3DOK4       14.64 4.116188 5 1.535491  9.1 19.4 11.9 15.8 17.0
## 3DOK5       23.98 3.777168 5 1.535491 17.7 27.9 24.3 24.8 25.2
## 3DOK7       19.92 1.130044 5 1.535491 18.6 21.0 18.8 20.5 20.7
## KOMPOSIT    18.70 1.601562 5 1.535491 16.9 20.8 17.3 19.1 19.4
## 
## $comparison
## NULL
## 
## $groups
##          Nitrogen groups
## 3DOK1       28.82      a
## 3DOK5       23.98      b
## 3DOK7       19.92     bc
## KOMPOSIT    18.70     cd
## 3DOK4       14.64     de
## 3DOK13      13.26      e
## 
## attr(,"class")
## [1] "group"

Interpretasi :

  1. 3DOK1 Paling tinggi, berbeda nyata dengan semua yang hurufnya bukan “a”.
  2. 3DOK2 Berbeda nyata dengan 3DOK1 (a) dan 3DOK13 (e), tetapi bisa sama dengan 3DOK7 (bc).
library(dplyr)
## Warning: package 'dplyr' was built under R version 4.3.3
## 
## Attaching package: 'dplyr'
## The following objects are masked from 'package:stats':
## 
##     filter, lag
## The following objects are masked from 'package:base':
## 
##     intersect, setdiff, setequal, union
library(ggplot2)
plot(hasil_Duncan, horiz = TRUE, las = 1)

Latihan Praktikum 1

Setelah memahami konsep dasar analisis ragam pada Rancangan Acak Lengkap (RAL), mahasiswa diminta untuk melakukan analisis uji lanjut pada data Praktikum 8 Sheet 1.

Data tersebut memuat dua variabel utama, yaitu:

  • Suhu sebagai faktor perlakuan
  • Kadar antosianin sebagai variabel respons

Latihan ini bertujuan agar mahasiswa mampu menentukan apakah terdapat perbedaan kadar antosianin pada setiap taraf suhu, serta mengidentifikasi pasangan perlakuan suhu yang berbeda nyata melalui uji lanjut.

Instruksi Latihan

  1. Impor data Praktikum 8 Sheet 1 ke dalam R.
  2. Pastikan variabel Suhu dibaca sebagai faktor perlakuan.
  3. Lakukan analisis ragam menggunakan model RAL.
  4. Jika hasil ANOVA menunjukkan pengaruh suhu yang signifikan, lanjutkan dengan uji lanjut.
  5. Gunakan salah satu metode uji lanjut, misalnya LSD, BNJ/Tukey, atau Duncan.
  6. Interpretasikan hasil uji lanjut berdasarkan kelompok perlakuan yang terbentuk.

Pertanyaan

Berdasarkan hasil analisis, apakah terdapat perbedaan kadar antosianin pada beberapa taraf suhu? Jika ada, suhu mana saja yang menghasilkan kadar antosianin yang berbeda nyata?

UJI LANJUT TERENCANA

CONTOH KONTRAS ORTHOGONAL

INPUT DATANYA

# Input Datanya
library(readxl)
## Warning: package 'readxl' was built under R version 4.3.3
data <- read_excel("Data_Kontras_Orthogonal.xlsx")

data$Pabrik <- as.factor(data$Pabrik)
data$Ulangan <- as.factor(data$Ulangan)
str(data)
## tibble [24 × 3] (S3: tbl_df/tbl/data.frame)
##  $ Pabrik : Factor w/ 6 levels "A","B","C","D",..: 1 1 1 1 2 2 2 2 3 3 ...
##  $ Ulangan: Factor w/ 4 levels "1","2","3","4": 1 2 3 4 1 2 3 4 1 2 ...
##  $ Bobot  : num [1:24] 6.78 6.51 6.92 7.08 6.88 6.13 7.21 7.84 9.8 9.72 ...
data
## # A tibble: 24 × 3
##    Pabrik Ulangan Bobot
##    <fct>  <fct>   <dbl>
##  1 A      1        6.78
##  2 A      2        6.51
##  3 A      3        6.92
##  4 A      4        7.08
##  5 B      1        6.88
##  6 B      2        6.13
##  7 B      3        7.21
##  8 B      4        7.84
##  9 C      1        9.8 
## 10 C      2        9.72
## # ℹ 14 more rows
#2. Buat Tabel Anova
model <- aov(Bobot ~ Pabrik, data = data)
summary(model)
##             Df Sum Sq Mean Sq F value   Pr(>F)    
## Pabrik       5 25.654   5.131   11.33 4.72e-05 ***
## Residuals   18  8.154   0.453                     
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Terbukti Pengaruh dan ingin dilakukan uji lanjut

#Install Terlebih Dahulu Pacakage (Emmeans)
library(emmeans)
## Welcome to emmeans.
## Caution: You lose important information if you filter this package's results.
## See '? untidy'
# 3. Rataan tiap pabrik
emm <- emmeans(model, ~ Pabrik)
emm
##  Pabrik emmean    SE df lower.CL upper.CL
##  A        6.82 0.337 18     6.12     7.53
##  B        7.01 0.337 18     6.31     7.72
##  C        9.80 0.337 18     9.09    10.50
##  D        7.50 0.337 18     6.80     8.21
##  E        6.88 0.337 18     6.17     7.58
##  F        7.96 0.337 18     7.25     8.67
## 
## Confidence level used: 0.95
# 4. Susun koefisien kontras ortogonal
kontras <- list(
  "H1: Dalam negeri vs luar negeri" = c( 2,  2, -1, -1, -1, -1),
  "H2: Pabrik A vs Pabrik B"        = c( 1, -1,  0,  0,  0,  0),
  "H3: C dan D vs E dan F"          = c( 0,  0,  1,  1, -1, -1),
  "H4: Pabrik C vs Pabrik D"        = c( 0,  0,  1, -1,  0,  0),
  "H5: Pabrik E vs Pabrik F"        = c( 0,  0,  0,  0,  1, -1)
)

# 5. Uji kontras
hasil_kontras <- contrast(emm, method = kontras, adjust = "none")
hasil_kontras
##  contrast                        estimate    SE df t.ratio p.value
##  H1: Dalam negeri vs luar negeri   -4.463 1.170 18  -3.828  0.0012
##  H2: Pabrik A vs Pabrik B          -0.192 0.476 18  -0.404  0.6906
##  H3: C dan D vs E dan F             2.462 0.673 18   3.659  0.0018
##  H4: Pabrik C vs Pabrik D           2.295 0.476 18   4.822  0.0001
##  H5: Pabrik E vs Pabrik F          -1.083 0.476 18  -2.274  0.0354
summary(hasil_kontras, infer = TRUE)
##  contrast                        estimate    SE df lower.CL upper.CL t.ratio
##  H1: Dalam negeri vs luar negeri   -4.463 1.170 18    -6.91  -2.0133  -3.828
##  H2: Pabrik A vs Pabrik B          -0.192 0.476 18    -1.19   0.8074  -0.404
##  H3: C dan D vs E dan F             2.462 0.673 18     1.05   3.8766   3.659
##  H4: Pabrik C vs Pabrik D           2.295 0.476 18     1.30   3.2949   4.822
##  H5: Pabrik E vs Pabrik F          -1.083 0.476 18    -2.08  -0.0826  -2.274
##  p.value
##   0.0012
##   0.6906
##   0.0018
##   0.0001
##   0.0354
## 
## Confidence level used: 0.95

CARA 2 AGAR HASIL SEPERTI TABEL ANOVA

# 1. Ambil rataan tiap perlakuan
mean_pabrik <- tapply(data$Bobot, data$Pabrik, mean)
mean_pabrik
##      A      B      C      D      E      F 
## 6.8225 7.0150 9.7975 7.5025 6.8775 7.9600
# 2. Ambil KT Galat dari ANOVA
anova_model <- anova(model)
KT_galat <- anova_model["Residuals", "Mean Sq"]
db_galat <- anova_model["Residuals", "Df"]

# 3. Banyak ulangan tiap perlakuan
r <- 4
# 4. Koefisien kontras ortogonal
kontras <- list(
  "H1: Dalam negeri vs luar negeri" = c( 2,  2, -1, -1, -1, -1),
  "H2: Pabrik A vs Pabrik B"        = c( 1, -1,  0,  0,  0,  0),
  "H3: C dan D vs E dan F"          = c( 0,  0,  1,  1, -1, -1),
  "H4: Pabrik C vs Pabrik D"        = c( 0,  0,  1, -1,  0,  0),
  "H5: Pabrik E vs Pabrik F"        = c( 0,  0,  0,  0,  1, -1)
)
# 5. Fungsi menghitung JK, KT, F hitung, dan p-value
hitung_kontras <- function(nama, c) {
  
  L <- sum(c * mean_pabrik)
  JK <- r * (L^2) / sum(c^2)
  db <- 1
  KT <- JK / db
  F_hitung <- KT / KT_galat
  p_value <- pf(F_hitung, db, db_galat, lower.tail = FALSE)
  
  data.frame(
    Sumber_Keragaman = nama,
    db = db,
    JK = JK,
    KT = KT,
    F_hitung = F_hitung,
    p_value = p_value
  )
}
# 6. Tabel ANOVA kontras ortogonal
tabel_kontras <- bind_rows(
  mapply(
    hitung_kontras,
    nama = names(kontras),
    c = kontras,
    SIMPLIFY = FALSE
  )
)

# 7. Rapikan angka
tabel_kontras <- tabel_kontras %>%
  mutate(
    JK = round(JK, 4),
    KT = round(KT, 4),
    F_hitung = round(F_hitung, 4),
    p_value = round(p_value, 4)
  )

tabel_kontras
##                  Sumber_Keragaman db      JK      KT F_hitung p_value
## 1 H1: Dalam negeri vs luar negeri  1  6.6380  6.6380  14.6526  0.0012
## 2        H2: Pabrik A vs Pabrik B  1  0.0741  0.0741   0.1636  0.6906
## 3          H3: C dan D vs E dan F  1  6.0639  6.0639  13.3854  0.0018
## 4        H4: Pabrik C vs Pabrik D  1 10.5340 10.5340  23.2528  0.0001
## 5        H5: Pabrik E vs Pabrik F  1  2.3436  2.3436   5.1733  0.0354
# 8. Tambahkan baris galat
tabel_galat <- data.frame(
  Sumber_Keragaman = "Galat",
  db = db_galat,
  JK = round(anova_model["Residuals", "Sum Sq"], 4),
  KT = round(KT_galat, 4),
  F_hitung = NA,
  p_value = NA
)

# 9. Gabungkan tabel kontras dan galat
tabel_anova_kontras <- bind_rows(tabel_kontras, tabel_galat)

tabel_anova_kontras
##                  Sumber_Keragaman db      JK      KT F_hitung p_value
## 1 H1: Dalam negeri vs luar negeri  1  6.6380  6.6380  14.6526  0.0012
## 2        H2: Pabrik A vs Pabrik B  1  0.0741  0.0741   0.1636  0.6906
## 3          H3: C dan D vs E dan F  1  6.0639  6.0639  13.3854  0.0018
## 4        H4: Pabrik C vs Pabrik D  1 10.5340 10.5340  23.2528  0.0001
## 5        H5: Pabrik E vs Pabrik F  1  2.3436  2.3436   5.1733  0.0354
## 6                           Galat 18  8.1544  0.4530       NA      NA

LATIHAN TIME

DATA ADA DI SHEET 2 FILE KONTRAS ORTHOGONAL

KONTRAS POLYNOMIAL

Koefisien pada tabel di atas merupakan koefisien polinomial orthogonal asli (standar) yang telah ditetapkan secara teoritis. Oleh karena itu, koefisien tersebut harus tetap digunakan apa adanya dan tidak diubah, selama jumlah perlakuan serta jarak antar taraf perlakuan sesuai.

Contoh Kontras Polynomial

library(readxl)
library(dplyr)
##install dulu knitr karena kaian belum pernah menggunakanya
library(knitr)

1. # Input Data
## [1] 1
# Baca data
data1 <- read_excel("Data_Kontras_Orthogonal.xlsx", sheet = "Polynomial1")
data1
## # A tibble: 30 × 3
##    Kelompok Jarak Produksi
##       <dbl> <dbl>    <dbl>
##  1        1    18     33.6
##  2        1    24     31.1
##  3        1    30     33  
##  4        1    36     28.4
##  5        1    42     31.4
##  6        2    18     37.1
##  7        2    24     34.5
##  8        2    30     29.5
##  9        2    36     29.9
## 10        2    42     28.3
## # ℹ 20 more rows
# Pastikan faktor
data1$Kelompok <- as.factor(data1$Kelompok)
data1$Jarak <- factor(data1$Jarak, levels = c(18, 24, 30, 36, 42))
# 2. Buatkan Tabel Anovanya
# ANOVA RAKL
model1 <- aov(Produksi ~ Kelompok + Jarak, data = data1)

# Tabel ANOVA
summary(model1)
##             Df Sum Sq Mean Sq F value   Pr(>F)    
## Kelompok     5   5.41   1.082   0.293 0.911327    
## Jarak        4 125.66  31.415   8.500 0.000354 ***
## Residuals   20  73.92   3.696                     
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
# 3. Lakukan Kontras Polynomial
# Rataan tiap jarak tanam
mean_jarak <- tapply(data1$Produksi, data1$Jarak, mean)
mean_jarak
##       18       24       30       36       42 
## 35.15000 31.63333 30.16667 29.53333 30.03333
# Ambil tabel ANOVA
anova_model <- anova(model1)

# Ambil KT Galat dan db Galat
KT_galat <- anova_model["Residuals", "Mean Sq"]
db_galat <- anova_model["Residuals", "Df"]

# Banyak kelompok/blok
r <- length(unique(data1$Kelompok))

# Koefisien polinomial ortogonal untuk 5 taraf perlakuan
kontras_poly <- list(
  "Linear"    = c(-2, -1,  0,  1,  2),
  "Kuadratik" = c( 2, -1, -2, -1,  2),
  "Kubik"     = c(-1,  2,  0, -2,  1),
  "Kuartik"   = c( 1, -4,  6, -4,  1)
)

# Fungsi menghitung JK, KT, F hitung, dan p-value
hitung_kontras <- function(nama, c) {
  
  L <- sum(c * mean_jarak)
  JK <- r * (L^2) / sum(c^2)
  db <- 1
  KT <- JK / db
  F_hitung <- KT / KT_galat
  p_value <- pf(F_hitung, db, db_galat, lower.tail = FALSE)
  
  data.frame(
    Sumber_Keragaman = nama,
    db = db,
    JK = JK,
    KT = KT,
    F_hitung = F_hitung,
    p_value = p_value
  )
}

# Membuat tabel kontras
tabel_kontras <- bind_rows(
  mapply(
    hitung_kontras,
    nama = names(kontras_poly),
    c = kontras_poly,
    SIMPLIFY = FALSE
  )
)

# Rapikan hasil
tabel_kontras <- tabel_kontras %>%
  mutate(
    JK = round(JK, 4),
    KT = round(KT, 4),
    F_hitung = round(F_hitung, 4),
    p_value = round(p_value, 4)
  )

tabel_kontras
##   Sumber_Keragaman db      JK      KT F_hitung p_value
## 1           Linear  1 91.2667 91.2667  24.6938  0.0001
## 2        Kuadratik  1 33.6933 33.6933   9.1163  0.0068
## 3            Kubik  1  0.5042  0.5042   0.1364  0.7158
## 4          Kuartik  1  0.1972  0.1972   0.0533  0.8197