Statistik deskriptif menggambar, uji statistik memutuskan. Pada bagian ini kita mengeksekusi rencana Bagian 1 langkah demi langkah: memeriksa asumsi normalitas lewat momen g1/g2, menguji H0 lewat ANOVA manual, lalu — karena H0 ditolak — mencari pasangan mana yang berbeda lewat uji t berpasangan dengan koreksi Bonferroni. Setiap uji dibangun dari rumus dasar dengan operasi inti MATLAB/Octave, termasuk p-value yang dihitung persis tanpa Statistics Toolbox.

Yang kita bangun pada bagian ini:

  • ujiNormalitas.m — skewness g1 dan kurtosis g2 dari rumus momen,
  • ujiAnova.m — tabel ANOVA lengkap: SS, df, MS, F, p-value, nilai kritis, eta-squared,
  • ujiPostHoc.m — uji t pooled untuk tiga pasangan dengan ambang Bonferroni.

Uji Normalitas lewat Momen #

Uji t dan ANOVA mengandaikan data berdistribusi mendekati normal. Untuk n = 10, uji formal seperti Shapiro-Wilk memang ideal, tetapi hidup di Statistics Toolbox. Kita memakai pemeriksaan ringan yang bisa dihitung manual: kemiringan (skewness) g1 dan keruncingan (kurtosis) g2 dari momen sampel:

  • m2 = mean((x-mean)^2), m3 = mean((x-mean)^3), m4 = mean((x-mean)^4)
  • g1 = m3 / m2^1.5 — distribusi normal simetris, g1 = 0
  • g2 = m4 / m2^2 - 3 — ekces kurtosis; normal memberi g2 = 0

Berkas src/ujiNormalitas.m:

matlab 
function N = ujiNormalitas(data, cetak)
% UJINORMALITAS Uji ringan normalitas lewat skewness dan kurtosis.
%   N = ujiNormalitas(data)        : hitung dan cetak tabel.
%   N = ujiNormalitas(data, false) : hitung tanpa mencetak.
%
%   Rumus momen klasik untuk satu vektor x dengan mean m:
%     m2 = mean((x-m)^2)   m3 = mean((x-m)^3)   m4 = mean((x-m)^4)
%     g1 = m3 / m2^1.5     (skewness)
%     g2 = m4 / m2^2 - 3   (excess kurtosis; normal = 0)
%
%   Kriteria ringan: |g1| < 0.5 dianggap cukup simetris untuk
%   melanjutkan ke ANOVA. Ini bukan uji formal - uji resmi seperti
%   Shapiro-Wilk atau Lilliefors (Statistics Toolbox MATLAB) lebih
%   tepat untuk sampel kecil; keterbatasan ini dibahas di Bagian 3.

if nargin < 2 || isempty(cetak)
    cetak = true;
end

k = numel(data.label);
N = repmat(struct('g1', 0, 'g2', 0, 'cukupNormal', false), 1, k);

if cetak
    fprintf('%-6s %9s %9s   %s\n', 'Metode', 'g1', 'g2', ...
        '|g1| < 0.5');
end

for i = 1:k
    x = data.waktu(i, :);
    m = mean(x);
    m2 = mean((x - m).^2);
    m3 = mean((x - m).^3);
    m4 = mean((x - m).^4);
    N(i).g1 = m3 / m2^1.5;
    N(i).g2 = m4 / m2^2 - 3;
    N(i).cukupNormal = abs(N(i).g1) < 0.5;
    if cetak
        if N(i).cukupNormal
            ket = 'cukup normal';
        else
            ket = 'kurang normal';
        end
        fprintf('%-6s %9.3f %9.3f   %s\n', data.label{i}, ...
            N(i).g1, N(i).g2, ket);
    end
end
end

Bedah singkatnya: momen m2-m4 dihitung dari deviasi x - m dengan operator .^ per elemen (baris 31-33); g1 adalah momen ketiga ternormalisasi pangkat 1,5 dari m2, g2 momen keempat ternormalisasi m2 kuadrat dikurangi 3 (baris 34-35); keputusan ringan abs(g1) < 0.5 disimpan sebagai logical di field cukupNormal (baris 36) sehingga orkestrator bisa menghentikan analisis bila ada grup yang melenceng. Menjalankannya:

matlab 
% driver: uji normalitas seluruh metode
data = buatDataEksperimen();
N = ujiNormalitas(data);

Output yang Dihasilkan:

text 
Metode         g1        g2   |g1| < 0.5
A          0.151    -1.190   cukup normal
B          0.029    -1.627   cukup normal
C          0.082    -1.440   cukup normal

Ketiga g1 kecil — sebaran tiap metode nyaris simetris, rencana lanjut ke ANOVA. Angka g2 yang serba negatif bukan tanda bahaya: pada sampel kecil, estimasi kurtosis memang terdorong turun — rata-rata g2 dari sampel normal berukuran n hanya sekitar -6/(n+1) = -0.55 untuk n = 10, dengan sebaran besar. Itulah alasan keputusan kita memakai g1 saja, dan g2 dilaporkan sebagai informasi bentuk.

ANOVA Satu Arah dari Rumus Dasar #

ANOVA membedah total variasi menjadi dua sumber, seperti chapter Analisis Variansi: variasi antara kelompok (mean metode menjauh dari mean keseluruhan) dan variasi dalam kelompok (replikat menjauh dari mean metodenya). Statistik F adalah rasio kedua rata-rata variasi itu — bila H0 benar, F berkisar di 1.

BesaranRumusNilai df
SSbetweenn * sum((meanGrup - grandMean)^2)k - 1 = 2
SSwithinsum((x - meanGrup)^2) untuk semua pengamatanN - k = 27
MSB / MSWSSbetween/dfB dan SSwithin/dfW-
FMSB / MSW-
eta^2SSbetween / (SSbetween + SSwithin)-

Berkas src/ujiAnova.m:

matlab 
function A = ujiAnova(data, alpha, cetak)
% UJIANOVA One-way ANOVA manual untuk kelompok seimbang.
%   A = ujiAnova(data)               : alpha 0.05, cetak tabel ANOVA.
%   A = ujiAnova(data, alpha)        : alpha kustom.
%   A = ujiAnova(data, alpha, false) : tanpa cetak.
%
%   Hipotesis: H0 mean semua kelompok sama, H1 minimal satu berbeda.
%   ANOVA membandingkan variabilitas ANTARA kelompok (SSbetween) dengan
%   variabilitas DALAM kelompok (SSwithin) lewat rasio F. Bila H0 benar,
%   F sekitar 1; F besar berarti perbedaan antar mean melampaui derau
%   replikasi.
%
%   p-value dan nilai kritis dihitung dari betainc/betaincinv yang ada
%   di inti MATLAB maupun Octave, sehingga tidak membutuhkan fcdf dari
%   Statistics Toolbox.

if nargin < 2 || isempty(alpha)
    alpha = 0.05;
end
if nargin < 3 || isempty(cetak)
    cetak = true;
end

X = data.waktu;
k = size(X, 1);              % jumlah kelompok
n = size(X, 2);              % replikat per kelompok (desain seimbang)
N = k * n;

mk = mean(X, 2);             % mean per kelompok (vektor kolom)
mg = mean(X(:));             % mean keseluruhan (grand mean)

% SS antara: sebaran mean kelompok terhadap mean keseluruhan.
A.ssBetween = n * sum((mk - mg).^2);
% SS dalam: sebaran tiap pengamatan terhadap mean kelompoknya.
% X - mk bekerja lewat implicit expansion (R2016b+).
deviasi = X - mk;
A.ssWithin = sum(deviasi(:).^2);

A.dfBetween = k - 1;
A.dfWithin = N - k;
A.msBetween = A.ssBetween / A.dfBetween;
A.msWithin = A.ssWithin / A.dfWithin;
A.F = A.msBetween / A.msWithin;

% p-value = P(F >= Fobs). CDF distribusi F lewat incomplete beta:
%   CDF_F(f; d1, d2) = betainc(d1*f / (d1*f + d2), d1/2, d2/2)
A.p = 1 - betainc(A.dfBetween * A.F / ...
    (A.dfBetween * A.F + A.dfWithin), A.dfBetween/2, A.dfWithin/2);

% Nilai kritis F: kebalikan CDF yang sama.
u = betaincinv(1 - alpha, A.dfBetween/2, A.dfWithin/2);
A.Fkritis = A.dfWithin * u / (A.dfBetween * (1 - u));

A.etaSquared = A.ssBetween / (A.ssBetween + A.ssWithin);
A.signifikan = A.p < alpha;

if cetak
    fprintf('Sumber        SS    df       MS       F        p\n');
    fprintf('Antara  %8.2f %5d %8.2f %7.2f %8.2e\n', ...
        A.ssBetween, A.dfBetween, A.msBetween, A.F, A.p);
    fprintf('Dalam   %8.2f %5d %8.2f\n', ...
        A.ssWithin, A.dfWithin, A.msWithin);
    fprintf('Total   %8.2f %5d\n', ...
        A.ssBetween + A.ssWithin, A.dfBetween + A.dfWithin);
    fprintf('F kritis (df %d, %d; alpha %.2f) = %.2f, eta^2 = %.3f\n', ...
        A.dfBetween, A.dfWithin, alpha, A.Fkritis, A.etaSquared);
    if A.signifikan
        fprintf('Keputusan: p < alpha -> H0 ditolak; minimal satu metode berbeda.\n');
    else
        fprintf('Keputusan: p >= alpha -> H0 tidak ditolak.\n');
    end
end
end

Mari kita bedah baris-baris kode di atas:

  1. Baris 29-30 (mk = mean(X, 2); mg = mean(X(:));): mean(X, 2) merata-ratakan per baris (per metode) menghasilkan vektor kolom 3x1; mean(X(:)) merata-ratakan seluruh 30 pengamatan menjadi grand mean.
  2. Baris 33 (A.ssBetween = n * sum((mk - mg).^2);): deviasi mean tiap metode terhadap grand mean dikuadratkan lalu dikalikan n — setiap mean grup mewakili 10 pengamatan.
  3. Baris 36-37 (deviasi = X - mk; A.ssWithin = sum(deviasi(:).^2);): X - mk mengurangkan mean metode dari tiap replikat lewat implicit expansion (matriks 3x10 dikurangi 3x1). Operator (:) meratakan deviasi menjadi kolom sebelum dijumlahkan — bentuk eksplisit yang setara sum(..., 'all').
  4. Baris 47-48 (A.p = 1 - betainc(...)): p-value = P(F sama besar atau lebih besar dari hasil amatan). CDF distribusi F dengan df (d1, d2) identik dengan incomplete beta function betainc(d1*f/(d1*f + d2), d1/2, d2/2) — fungsi inti yang hadir di MATLAB maupun Octave.
  5. Baris 51-52 (u = betaincinv(...); A.Fkritis = ...): nilai kritis adalah kebalikan CDF pada 1 - alpha. Hasil transformasinya kembali ke skala F lewat dfW*u/(dfB*(1-u)).
  6. Baris 54-55 (A.etaSquared ...; A.signifikan ...): dua ringkasan keputusan — effect size dan logical signifikansi — siap dipakai laporan tanpa menghitung ulang.

Menjalankannya:

matlab 
% driver: ANOVA satu arah
data = buatDataEksperimen();
A = ujiAnova(data);

Output yang Dihasilkan:

text 
Sumber        SS    df       MS       F        p
Antara    260.00     2   130.00   25.30 6.46e-07
Dalam     138.74    27     5.14
Total     398.74    29
F kritis (df 2, 27; alpha 0.05) = 3.35, eta^2 = 0.652
Keputusan: p < alpha -> H0 ditolak; minimal satu metode berbeda.

p-value Tanpa Toolbox #

Chapter Analisis Variansi berhenti pada perbandingan F dengan nilai kritis tabel karena fcdf hidup di Statistics Toolbox (dan tidak ada di Octave — exist('fcdf') mengembalikan 0). Proyek ini melangkah lebih jauh: karena betainc/betaincinv tersedia di inti kedua mesin, p-value dihitung persis. Benarkanya mudah dicek terhadap angka tabel klasik:

matlab 
% verifikasi identitas betainc terhadap nilai tabel klasik
fprintf('p di F = 3.35 (df 2, 27) : %.4f\n', ...
    1 - betainc(2*3.35/(2*3.35 + 27), 1, 13.5));
u = betaincinv(0.95, 1, 13.5);
fprintf('F kritis dari betaincinv : %.4f\n', 27*u/(2*(1 - u)));
fprintf('p dua-ekor t = 2.101 (df 18): %.4f\n', ...
    betainc(18/(18 + 2.101^2), 9, 0.5));

Output yang Dihasilkan:

text 
p di F = 3.35 (df 2, 27) : 0.0502
F kritis dari betaincinv : 3.3541
p dua-ekor t = 2.101 (df 18): 0.0500

Nilai tabel klasik F(2,27) di alpha 0.05 adalah 3.35 — betaincinv memberi 3.3541; memasukkan 3.35 ke identitas CDF mengembalikan p = 0.0502, nyaris tepat alpha. Identitas t juga cocok: t = 2.101 pada df 18 adalah titik kritis dua-ekor 0.05, dan betainc menghitung p = 0.0500. Rumus-rumus itu sama yang dipakai fcdf/tcdf di balik layar.

Membaca Hasil ANOVA #

F = 25.30 jauh melampaui nilai kritis 3.35, dan p = 6.46e-07 jauh di bawah alpha: H0 ditolak dengan bukti sangat kuat. Variasi antara metode (MSB = 130.00) adalah 25 kali variasi derau dalam metode (MSW = 5.14). Effect size menegaskan besarnya: eta^2 = 0.652 berarti 65.2% variasi waktu produksi dijelaskan oleh pilihan metode — jauh di atas patokan Cohen untuk efek besar:

eta^2Interpretasi
0.01kecil
0.06sedang
0.14besar
0.652 (proyek ini)sangat besar

Posthoc dengan Koreksi Bonferroni #

ANOVA yang signifikan hanya berkata “ada perbedaan” — pasangan mana, tidak. Pertanyaan itu dijawab uji t untuk tiap pasangan, tetapi tiap uji membawa peluang kesalahan tipe I sendiri, dan peluang setidaknya satu kesalahan (family-wise error, FWE) menggelembung:

matlab 
% family-wise error tanpa koreksi
for m = [1 2 3 10]
    fprintf('m = %2d perbandingan: FWE = %.3f\n', m, 1 - (1 - 0.05)^m);
end

Output yang Dihasilkan:

text 
m =  1 perbandingan: FWE = 0.050
m =  2 perbandingan: FWE = 0.098
m =  3 perbandingan: FWE = 0.143
m = 10 perbandingan: FWE = 0.401

Tiga perbandingan tanpa koreksi menaikkan peluang temuan palsu menjadi 14.3%. Bonferroni menutupnya dengan cara paling sederhana: bagi alpha dengan banyaknya perbandingan — 0.05/3 = 0.0167. Berkas src/ujiPostHoc.m:

matlab 
function P = ujiPostHoc(data, alpha, cetak)
% UJIPOSTHOC Uji t berpasangan dengan koreksi Bonferroni.
%   P = ujiPostHoc(data)              : alpha 0.05, cetak tabel.
%   P = ujiPostHoc(data, alpha)       : alpha kustom.
%   P = ujiPostHoc(data, alpha, false): tanpa cetak.
%
%   Setelah ANOVA signifikan, kita cari pasangan mana yang berbeda.
%   k kelompok menghasilkan k*(k-1)/2 perbandingan; bila setiap uji t
%   memakai alpha sendiri, peluang minimal satu kesalahan tipe I
%   (family-wise error) menggelembung: 1-(1-alpha)^m. Bonferroni
%   menutupnya dengan ambang alpha/m.
%
%   Uji t dua sampel memakai varians gabungan (pooled):
%     sp2 = ((n1-1)s1^2 + (n2-1)s2^2) / (n1+n2-2)
%     t   = (mean1 - mean2) / sqrt(sp2 * (1/n1 + 1/n2))
%   p dua-ekor dihitung dari betainc:
%     p = P(|T| >= |t|) = betainc(df/(df+t^2), df/2, 1/2)

if nargin < 2 || isempty(alpha)
    alpha = 0.05;
end
if nargin < 3 || isempty(cetak)
    cetak = true;
end

X = data.waktu;
k = size(X, 1);
n = size(X, 2);

% Daftar semua pasangan (i, j) dengan i < j.
pasangan = zeros(k * (k - 1) / 2, 2);
baris = 0;
for i = 1:k
    for j = i + 1:k
        baris = baris + 1;
        pasangan(baris, :) = [i, j];
    end
end
m = size(pasangan, 1);            % jumlah perbandingan
alphaDikoreksi = alpha / m;

P = repmat(struct('pasangan', '', 'selisih', 0, 't', 0, 'df', 0, ...
    'p', 1, 'signifikan', false), 1, m);

for c = 1:m
    i = pasangan(c, 1);
    j = pasangan(c, 2);
    xa = X(i, :);
    xb = X(j, :);
    selisih = mean(xa) - mean(xb);
    sp2 = ((n - 1) * var(xa) + (n - 1) * var(xb)) / (2 * n - 2);
    se = sqrt(sp2 * (1 / n + 1 / n));
    df = 2 * n - 2;
    tc = selisih / se;
    pc = betainc(df / (df + tc^2), df / 2, 0.5);

    P(c).pasangan = sprintf('%s - %s', data.label{i}, data.label{j});
    P(c).selisih = selisih;
    P(c).t = tc;
    P(c).df = df;
    P(c).p = pc;
    P(c).signifikan = pc < alphaDikoreksi;
end

if cetak
    fprintf('Jumlah perbandingan: %d, alpha Bonferroni = %.2f / %d = %.4f\n', ...
        m, alpha, m, alphaDikoreksi);
    fprintf('%-8s %8s %10s %10s   %s\n', 'Pasangan', 'selisih', ...
        't', 'p', 'keputusan');
    for c = 1:m
        if P(c).signifikan
            ket = 'signifikan';
        else
            ket = 'tidak signifikan';
        end
        fprintf('%-8s %+8.2f %+10.2f %10.4g   %s\n', P(c).pasangan, ...
            P(c).selisih, P(c).t, P(c).p, ket);
    end
end
end

Mari kita bedah baris-baris kode di atas:

  1. Baris 31-38 (loop ganda for i ... for j = i+1 ...): daftar pasangan dibangun generik — k grup menghasilkan k(k-1)/2 baris; untuk k = 3 pasangannya (1,2), (1,3), (2,3). Tidak ada pasangan terbalik atau pasangan diri-sendiri karena j selalu mulai dari i+1.
  2. Baris 39-40 (m = size(pasangan, 1); alphaDikoreksi = alpha / m;): inti Bonferroni — ambang dibagi jumlah perbandingan, bukan angka ajaib yang di-hardcode.
  3. Baris 51 (sp2 = ((n - 1) * var(xa) + ...): varians gabungan (pooled) membobot SS kedua grup; dengan desain seimbang sederhananya menjadi rata-rata dua varians.
  4. Baris 52-54 (se, df, tc): standard error selisih mean memakai sp2 dan kedua n; df = n1 + n2 - 2 = 18; t = selisih dibagi se.
  5. Baris 55 (pc = betainc(df / (df + tc^2), df / 2, 0.5);): p dua-ekor dari identitas t-beta — sama persis dengan yang diverifikasi di blok sebelumnya (t = 2.101 memberi p = 0.0500).
  6. Baris 62 (P(c).signifikan = pc < alphaDikoreksi;): keputusan membandingkan p mentah dengan ambang terkoreksi — ekuivalen dengan membandingkan p*m terhadap alpha.

Menjalankannya:

matlab 
% driver: posthoc Bonferroni
data = buatDataEksperimen();
P = ujiPostHoc(data);

Output yang Dihasilkan:

text 
Jumlah perbandingan: 3, alpha Bonferroni = 0.05 / 3 = 0.0167
Pasangan  selisih          t          p   keputusan
A - B       +5.00      +5.80  1.698e-05   signifikan
A - C       -2.00      -1.86    0.07994   tidak signifikan
B - C       -7.00      -6.45  4.565e-06   signifikan

Tanda selisih mengikuti urutan pasangan: A - B = +5 berarti A lima menit lebih lambat dari B. Dua pasangan yang melibatkan B signifikan bahkan setelah ambang Bonferroni mengetat — p keduanya berorde 1e-05, ribuan kali lebih kecil dari 0.0167. Pasangan A - C tidak signifikan: p = 0.07994 gagal melampaui alpha 0.05 biasa, apalagi ambang terkoreksi.

Perbandingan dengan Statistics Toolbox #

Pemilik MATLAB dengan Statistics and Machine Learning Toolbox bisa memverifikasi seluruh pipeline dengan tiga fungsi:

matlab 
% MATLAB dengan Statistics and Machine Learning Toolbox
[pTabel, tbl, stats] = anova1(waktu30, grup30, 'off');
c = multcompare(stats, 'CType', 'bonferroni');
[h, pAB] = ttest2(xA, xB);

anova1 menerima data vektor panjang plus vektor grup, mengembalikan p-value, tabel ANOVA, dan struct stats; multcompare dengan 'CType', 'bonferroni' melakukan uji banding ber-koreksi Bonferroni; ttest2 menguji satu pasangan dan mengembalikan keputusan h serta p-value p. Nilai F, df, dan p yang dihasilkan bertumpang tindih dengan implementasi manual kita — tabel ANOVA ujiAnova dirancang mengikuti kolom-kolom tbl milik anova1.

Tiga angka kini menggantung menunggu penafsiran: B tercepat, A dan C tidak terpisahkan. Bagian 4 merangkai semuanya menjadi kesimpulan yang didukung data, laporan tertulis, dan daftar keterbatasan yang jujur.