library(MASS)

1 Perbandingan Dua Vektor Nilai Tengah Berpasangan

1.1 Input Data

Suatu studi akan membandingkan dua jenis pupuk yang akan digunakan pada suatu lahan tanah yang sama. Kemudian, akan diukur pengaruhnya pada tanaman (Pupuk A = kontrol lama, Pupuk B = formulasi baru). Data diasumsikan multivariat normal, kedua sampel berpasangan karena akan diujikan pada lahan tanah yang sama, dan kedua populasi memiliki kovarians yang sama. Untuk setiap tanaman diukur dua metrik setelah 8 minggu:

Height = tinggi tumbuhan (cm)

Chloro = indeks klorofil (skala 0–100)

Hipotesis:

\(H_0: \mu_A = \mu_B\)

\(H_1: \mu_A ≠ \mu_B\)

dengan \(\mu = [Height, Chloro]^T\)

set.seed(80)

n1 <- n2 <- 30
mu_A <- c(55,48)
mu_B <- c(58,51)
sigma <- matrix(c(25,4,
                  4,9), nrow=2, byrow=TRUE)

X <- mvrnorm(n=n1, mu = mu_A, Sigma = sigma)
Y <- mvrnorm(n=n2, mu = mu_B, Sigma = sigma)

df <- data.frame(
  Pupuk = rep(c("A","B"), each=30),
  Height = c(X[,1],Y[,1]),
  Chloro = c(X[,2],Y[,2])
)

df
##    Pupuk   Height   Chloro
## 1      A 54.00692 55.57549
## 2      A 51.28501 46.19415
## 3      A 54.31736 45.01633
## 4      A 52.06249 50.48225
## 5      A 53.92232 47.38684
## 6      A 57.18095 46.78990
## 7      A 54.00572 51.99143
## 8      A 55.13819 52.32783
## 9      A 58.24816 44.95969
## 10     A 53.56519 52.48927
## 11     A 56.02405 46.26567
## 12     A 60.54250 51.22355
## 13     A 61.25595 53.29287
## 14     A 48.27271 47.41084
## 15     A 55.16066 47.11228
## 16     A 51.60787 48.51845
## 17     A 56.46204 48.30581
## 18     A 55.87080 47.70652
## 19     A 48.26533 48.97215
## 20     A 57.46212 48.52388
## 21     A 59.79895 49.49578
## 22     A 57.91177 51.03246
## 23     A 51.63520 45.79091
## 24     A 57.49883 45.75733
## 25     A 51.19101 49.56069
## 26     A 57.00634 45.38665
## 27     A 43.09978 47.24048
## 28     A 41.51653 48.86308
## 29     A 48.15568 43.50933
## 30     A 62.60775 48.72135
## 31     B 60.48178 51.81504
## 32     B 49.87814 52.77584
## 33     B 58.34507 58.07347
## 34     B 53.01056 51.60254
## 35     B 54.47587 49.84549
## 36     B 57.48962 51.50391
## 37     B 50.06472 44.69551
## 38     B 59.31940 50.75080
## 39     B 64.39305 49.66756
## 40     B 55.49234 50.78152
## 41     B 61.75740 44.13282
## 42     B 53.96986 47.59978
## 43     B 59.17478 50.94463
## 44     B 60.13357 55.29240
## 45     B 62.48521 50.64751
## 46     B 52.71049 48.49593
## 47     B 43.80085 51.53002
## 48     B 42.54269 43.03497
## 49     B 60.40866 53.10029
## 50     B 61.85072 53.97341
## 51     B 60.03753 52.41390
## 52     B 58.32000 51.97523
## 53     B 56.58864 49.06369
## 54     B 56.96013 51.81538
## 55     B 69.61117 51.17924
## 56     B 57.62131 51.81251
## 57     B 56.36232 50.42696
## 58     B 71.71458 59.11368
## 59     B 56.28708 52.27402
## 60     B 51.47664 52.25647

1.1.1 Pemisahan Kelompok Pupuk

colnames(X) <- colnames(Y) <- c("Height","Chloro")
X; Y
##         Height   Chloro
##  [1,] 54.00692 55.57549
##  [2,] 51.28501 46.19415
##  [3,] 54.31736 45.01633
##  [4,] 52.06249 50.48225
##  [5,] 53.92232 47.38684
##  [6,] 57.18095 46.78990
##  [7,] 54.00572 51.99143
##  [8,] 55.13819 52.32783
##  [9,] 58.24816 44.95969
## [10,] 53.56519 52.48927
## [11,] 56.02405 46.26567
## [12,] 60.54250 51.22355
## [13,] 61.25595 53.29287
## [14,] 48.27271 47.41084
## [15,] 55.16066 47.11228
## [16,] 51.60787 48.51845
## [17,] 56.46204 48.30581
## [18,] 55.87080 47.70652
## [19,] 48.26533 48.97215
## [20,] 57.46212 48.52388
## [21,] 59.79895 49.49578
## [22,] 57.91177 51.03246
## [23,] 51.63520 45.79091
## [24,] 57.49883 45.75733
## [25,] 51.19101 49.56069
## [26,] 57.00634 45.38665
## [27,] 43.09978 47.24048
## [28,] 41.51653 48.86308
## [29,] 48.15568 43.50933
## [30,] 62.60775 48.72135
##         Height   Chloro
##  [1,] 60.48178 51.81504
##  [2,] 49.87814 52.77584
##  [3,] 58.34507 58.07347
##  [4,] 53.01056 51.60254
##  [5,] 54.47587 49.84549
##  [6,] 57.48962 51.50391
##  [7,] 50.06472 44.69551
##  [8,] 59.31940 50.75080
##  [9,] 64.39305 49.66756
## [10,] 55.49234 50.78152
## [11,] 61.75740 44.13282
## [12,] 53.96986 47.59978
## [13,] 59.17478 50.94463
## [14,] 60.13357 55.29240
## [15,] 62.48521 50.64751
## [16,] 52.71049 48.49593
## [17,] 43.80085 51.53002
## [18,] 42.54269 43.03497
## [19,] 60.40866 53.10029
## [20,] 61.85072 53.97341
## [21,] 60.03753 52.41390
## [22,] 58.32000 51.97523
## [23,] 56.58864 49.06369
## [24,] 56.96013 51.81538
## [25,] 69.61117 51.17924
## [26,] 57.62131 51.81251
## [27,] 56.36232 50.42696
## [28,] 71.71458 59.11368
## [29,] 56.28708 52.27402
## [30,] 51.47664 52.25647

1.2 Menghitung Rataan dan Matriks Kovarians

1.2.1 Rataan Per Kelompok Pupuk

xbar <- colMeans(X)
ybar <- colMeans(Y)
S1 <- cov(X)
S2 <- cov(Y)

xbar; ybar
##   Height   Chloro 
## 54.16927 48.53011
##   Height   Chloro 
## 57.22547 51.08648

1.2.2 Kovarians Per Kelompok Pupuk

S1; S2
##           Height   Chloro
## Height 24.154179 2.481001
## Chloro  2.481001 7.931533
##           Height    Chloro
## Height 38.889709  9.794112
## Chloro  9.794112 11.556721

1.2.3 Kovarians Gabungan

Sp <- ((n1 - 1) * S1 + (n2 - 1) * S2) / (n1 + n2 - 2)
Sp
##           Height   Chloro
## Height 31.521944 6.137556
## Chloro  6.137556 9.744127

1.3 Uji T2 Hotelling Dua Populasi Sampel Berpasangan

diff <- xbar - ybar
T2 <- (n1 * n2) / (n1 + n2) * as.numeric(t(diff) %*% solve(Sp) %*% diff)
T2
## [1] 11.19408

1.3.1 Transformasi ke F-Stat

p = 2
Fstat <- ( (n1 + n2 - p - 1) / ( (n1 + n2 - 2) * p ) ) * T2
df1 <- p
df2 <- n1 + n2 - p - 1
p_value <- 1 - pf(Fstat, df1, df2)
Fstat; p_value
## [1] 5.500537
## [1] 0.006542351

1.4 Interpretasi

Hasil Uji T2 Hotelling yang dikonversi menjadi nilai F menunjukkan bahwa nilai p-value < 0.05. Maka tolak H0 pada taraf nyata 5%, cukup bukti untuk menyatakan terdapat perbedaan signifikan secara multivariat (Height dan Chloro) antara pupuk A dan pupuk B. Terlihat bahwa pupuk B menghasilkan tanaman yang lebih tinggi dengan kandungan klorofil yang lebih tinggi jika dibandingkan dengan tanaman yang diberi pupuk A.

1.5 Selang Kepercayaan

1.5.1 Simultan

C <- (1/n1 + 1/n2) * Sp
se <- sqrt(diag(C))
alpha <- 0.05

crit_simul <- sqrt(p*(n1+n2-1)/(n1+n2-p) * qf(1-alpha, p, n1+n2-p))
CI_simul <- cbind(
  lower = diff - crit_simul * se,
  upper = diff + crit_simul * se
)

CI_simul
##            lower      upper
## Height -6.729460  0.6170604
## Chloro -4.598663 -0.5140898

Pada taraf nyata 5%, selang kepercayaan simultan untuk height berada dalam rentang -6.729 sampai 0.617. Sedangkan selang kepercayaan simultan untuk chloro berada dalam rentang -4.598 sampai -0.514. Lewat selang kepercayaan simultan ini, dapat disimpulkan bahwa untuk variabel height, tidak terdapat perbedaan signifikan pada penggunaan pupuk A dan B karena memuat nilai 0 dalam selang tersebut. Sedangkan terdapat perbedaan yang signifikan pada variabel chloro antara penggunaan pupuk A dan B.

1.5.2 Bonferroni

crit_bonf <- qt(1-alpha/(2*p), df=n1+n2-2)
CI_bonf <- cbind(
  lower = diff - crit_bonf * se,
  upper = diff + crit_bonf * se
)

CI_bonf
##            lower      upper
## Height -6.391948  0.2795484
## Chloro -4.411010 -0.7017422

Pada taraf nyata 5%, selang kepercayaan Bonferroni untuk height berada dalam rentang -6.391 sampai 0.279. Sedangkan selang kepercayaan bonferroni untuk chloro berada dalam rentang -4.411 sampai -0.701. Lewat selang kepercayaan simultan ini, dapat disimpulkan bahwa untuk variabel height, tidak terdapat perbedaan signifikan pada penggunaan pupuk A dan B karena memuat nilai 0 dalam selang tersebut. Sedangkan terdapat perbedaan yang signifikan pada variabel chloro antara penggunaan pupuk A dan B.

2 Perbandingan Dua Vektor Nilai Tengah Saling Bebas

2.1 Input Data

Sebuah perusahaan teknologi ingin membandingkan efektivitas dua program pelatihan (Program X dan Program Y) untuk karyawan baru. Pada program pelatihan ini, diasumsikan data menyebar multivariate normal dan saling bebas (independen) karena akan mmasing-masing pelatihan akan diterapkan pada karyawan baru yang berbeda. Dua indikator kinerja diukur setelah pelatihan:

Speed = jumlah tugas yang dapat diselesaikan per jam.

Accuracy = persentase penyelesaian tanpa kesalahan.

\(H_0: \mu_X = \mu_Y\)

\(H_1: \mu_X ≠ \mu_Y\)

dengan \(\mu = [Speed, Accuracy]^T\)

set.seed(80)

n1 <- n2 <- 30
mu_X <- c(50, 85)
mu_Y <- c(55, 80)
sigma2 <- matrix(c(20, 5,
                   5, 15), nrow=2)

x <- mvrnorm(n1, mu_X, Sigma=sigma2)
y <- mvrnorm(n2, mu_Y, Sigma=sigma2)

colnames(x) <- colnames(y) <- c("Speed", "Accuracy")
df1 <- data.frame(
  group = factor(c(rep("Program_X", n1), rep("Program_Y", n2))),
  rbind(x, y)
)

df1
##        group    Speed Accuracy
## 1  Program_X 45.76225 93.24567
## 2  Program_X 47.34340 82.06595
## 3  Program_X 50.67272 81.48923
## 4  Program_X 46.18765 87.06171
## 5  Program_X 49.26844 84.04880
## 6  Program_X 52.55344 84.17845
## 7  Program_X 47.33010 89.22910
## 8  Program_X 48.23366 89.88360
## 9  Program_X 54.34490 82.38906
## 10 Program_X 46.70340 89.67902
## 11 Program_X 51.70944 83.30749
## 12 Program_X 53.73172 89.97057
## 13 Program_X 53.48785 92.46429
## 14 Program_X 44.01567 82.69115
## 15 Program_X 50.53769 84.04459
## 16 Program_X 46.62549 84.74967
## 17 Program_X 51.22275 85.70099
## 18 Program_X 50.93649 84.88453
## 19 Program_X 43.32534 84.43894
## 20 Program_X 52.05526 86.19044
## 21 Program_X 53.79814 87.85222
## 22 Program_X 51.37432 89.11174
## 23 Program_X 47.84486 81.69990
## 24 Program_X 53.30042 83.09926
## 25 Program_X 45.78243 85.81544
## 26 Program_X 53.00571 82.56318
## 27 Program_X 39.29030 81.23254
## 28 Program_X 37.11089 82.66280
## 29 Program_X 45.61502 78.29046
## 30 Program_X 56.74343 87.67275
## 31 Program_Y 56.94604 81.52153
## 32 Program_Y 46.68634 79.99960
## 33 Program_Y 52.22369 88.01106
## 34 Program_Y 50.10654 79.45246
## 35 Program_Y 52.23537 77.84261
## 36 Program_Y 54.30582 80.43960
## 37 Program_Y 50.39675 70.99054
## 38 Program_Y 56.33336 80.04409
## 39 Program_Y 61.51539 80.07360
## 40 Program_Y 52.76879 79.14062
## 41 Program_Y 61.49268 73.22550
## 42 Program_Y 52.74894 75.20208
## 43 Program_Y 56.11432 80.22585
## 44 Program_Y 55.10068 85.33291
## 45 Program_Y 59.31612 80.70418
## 46 Program_Y 51.18806 75.89767
## 47 Program_Y 41.59261 77.11419
## 48 Program_Y 44.14400 67.28636
## 49 Program_Y 56.31556 82.94386
## 50 Program_Y 57.27142 84.27567
## 51 Program_Y 56.27166 82.08374
## 52 Program_Y 54.87001 81.17126
## 53 Program_Y 54.53805 77.48430
## 54 Program_Y 53.67816 80.65887
## 55 Program_Y 65.69552 83.04637
## 56 Program_Y 54.29293 80.81769
## 57 Program_Y 53.73125 78.95652
## 58 Program_Y 64.17386 92.45314
## 59 Program_Y 52.85286 81.00788
## 60 Program_Y 48.39695 79.80934

2.1.1 Pemisahan Kelompok Pelatihan

x; y
##          Speed Accuracy
##  [1,] 45.76225 93.24567
##  [2,] 47.34340 82.06595
##  [3,] 50.67272 81.48923
##  [4,] 46.18765 87.06171
##  [5,] 49.26844 84.04880
##  [6,] 52.55344 84.17845
##  [7,] 47.33010 89.22910
##  [8,] 48.23366 89.88360
##  [9,] 54.34490 82.38906
## [10,] 46.70340 89.67902
## [11,] 51.70944 83.30749
## [12,] 53.73172 89.97057
## [13,] 53.48785 92.46429
## [14,] 44.01567 82.69115
## [15,] 50.53769 84.04459
## [16,] 46.62549 84.74967
## [17,] 51.22275 85.70099
## [18,] 50.93649 84.88453
## [19,] 43.32534 84.43894
## [20,] 52.05526 86.19044
## [21,] 53.79814 87.85222
## [22,] 51.37432 89.11174
## [23,] 47.84486 81.69990
## [24,] 53.30042 83.09926
## [25,] 45.78243 85.81544
## [26,] 53.00571 82.56318
## [27,] 39.29030 81.23254
## [28,] 37.11089 82.66280
## [29,] 45.61502 78.29046
## [30,] 56.74343 87.67275
##          Speed Accuracy
##  [1,] 56.94604 81.52153
##  [2,] 46.68634 79.99960
##  [3,] 52.22369 88.01106
##  [4,] 50.10654 79.45246
##  [5,] 52.23537 77.84261
##  [6,] 54.30582 80.43960
##  [7,] 50.39675 70.99054
##  [8,] 56.33336 80.04409
##  [9,] 61.51539 80.07360
## [10,] 52.76879 79.14062
## [11,] 61.49268 73.22550
## [12,] 52.74894 75.20208
## [13,] 56.11432 80.22585
## [14,] 55.10068 85.33291
## [15,] 59.31612 80.70418
## [16,] 51.18806 75.89767
## [17,] 41.59261 77.11419
## [18,] 44.14400 67.28636
## [19,] 56.31556 82.94386
## [20,] 57.27142 84.27567
## [21,] 56.27166 82.08374
## [22,] 54.87001 81.17126
## [23,] 54.53805 77.48430
## [24,] 53.67816 80.65887
## [25,] 65.69552 83.04637
## [26,] 54.29293 80.81769
## [27,] 53.73125 78.95652
## [28,] 64.17386 92.45314
## [29,] 52.85286 81.00788
## [30,] 48.39695 79.80934

2.2 Menghitung Rataan dan Matriks Kovarians

2.2.1 Rataan Per Kelompok Pelatihan

xbar1 <- colMeans(x)
ybar1 <- colMeans(y)
S1_1 <- cov(x)
S2_1 <- cov(y)

xbar1; ybar1
##    Speed Accuracy 
## 48.99711 85.39045
##    Speed Accuracy 
## 54.24346 79.90710

2.2.2 Kovarians Per Kelompok Pelatihan

S1_1; S2_1
##              Speed  Accuracy
## Speed    20.301029  3.915314
## Accuracy  3.915314 12.773196
##             Speed Accuracy
## Speed    27.74179 12.30737
## Accuracy 12.30737 22.22711

2.3 Kasus Ragam Sama

2.3.1 Matriks Kovarians Gabungan

Sp1 <- ((n1 - 1) * S1_1 + (n2 - 1) * S2_1) / (n1 + n2 - 2)
Sp1
##              Speed  Accuracy
## Speed    24.021411  8.111342
## Accuracy  8.111342 17.500154

2.3.2 Uji T2 Hotelling Dua Populasi Saling Bebas Ragam Sama

library(Hotelling)
## Loading required package: corpcor
t2 <- hotelling.test(x, y, var.equal=TRUE)
t2
## Test stat:  70.672 
## Numerator df:  2 
## Denominator df:  57 
## P-value:  1.372e-10

2.3.3 Interpretasi

Hasil Uji T2 Hotelling menunjukkan bahwa nilai p-value < 0.05. Maka tolak H0 pada taraf nyata 5%, cukup bukti untuk menyatakan terdapat perbedaan signifikan secara multivariat (Speed dan Accuracy) antara pelatihan X dan pelatihan Y. Terlihat bahwa pelatihan X memberikan nilai akurasi yang lebih tinggi, sedangkan pelatihan Y memberikan nilai speed / kecepatan yang lebih tinggi.

2.3.4 Selang Kepercayaan

2.3.4.1 Simultan

T.ci = function(mu_X, mu_Y, Sp1, n1, n2, avec=rep(1,length(mu)), level=0.95){
p = length(mu_X)
mu = mu_X-mu_Y
cval = qf(level, p, n1+n2-p-1) * p * (n1+n2-2) / (n1+n2-p-1)
zhat = crossprod(avec, mu)
zvar = crossprod(avec, Sp1 %*% avec)* (1/n1+1/n2)
const = sqrt(cval * zvar)
c(lower = zhat - const, upper = zhat + const)
}

T.ci(xbar1, ybar1, Sp1, n1,n2, avec=c(1,0),level=0.95) #Speed
##     lower     upper 
## -8.454903 -2.037797
T.ci(xbar1, ybar1, Sp1, n1,n2, avec=c(0,1),level=0.95) #Accuracy
##    lower    upper 
## 2.744735 8.221964

Pada taraf nyata 5%, selang kepercayaan simultan untuk speed berada dalam rentang -8.454 sampai -2.037. Sedangkan selang kepercayaan simultan untuk accuracy berada dalam rentang 2.744 sampai 8.221. Lewat selang kepercayaan simultan ini, dapat disimpulkan bahwa untuk variabel speed maupun accuracy, terdapat perbedaan signifikan antara pelatihan X dengan pelatihan Y.

2.3.4.2 Bonferroni

bon = function(mu_X, mu_Y , Sp1, n1, n2, alpha, k){
 p = length(mu_X)
 mu = mu_X-mu_Y
 lower = mu[k] - sqrt((Sp1[k,k]) *(1/n1+1/n2))* abs(qt(alpha/(2*p), df=n1+n2-2))
 upper = mu[k] + sqrt((Sp1[k,k]) *(1/n1+1/n2))* abs(qt(alpha/(2*p), df=n1+n2-2))
 ci = c(lower = lower,upper = upper)
 names(ci)= c("lower","upper")
 ci
}

bon(xbar1, ybar1, Sp1, n1, n2, 0.05, 1) #Speed
##     lower     upper 
## -8.158314 -2.334386
bon(xbar1, ybar1, Sp1, n1, n2, 0.05, 2) #Accuracy
##    lower    upper 
## 2.997884 7.968815

Pada taraf nyata 5%, selang kepercayaan simultan untuk speed berada dalam rentang -8.158 sampai -2.334. Sedangkan selang kepercayaan simultan untuk accuracy berada dalam rentang 2.997 sampai 7.968. Lewat selang kepercayaan bonferroni ini, dapat disimpulkan bahwa untuk variabel speed maupun accuracy, terdapat perbedaan signifikan antara pelatihan X dengan pelatihan Y.

2.4 Kasus Ragam Tidak Sama

2.4.1 Uji T2 Hotelling Dua Populasi Saling Bebas Ragam Tidak Sama

library(Hotelling)
t2_1 <- hotelling.test(x, y, var.equal=FALSE)
t2_1
## Test stat:  70.672 
## Numerator df:  2 
## Denominator df:  55.1267125271191 
## P-value:  1.738e-10

2.4.2 Interpretasi

Hasil Uji T2 Hotelling menunjukkan bahwa nilai p-value < 0.05. Maka tolak H0 pada taraf nyata 5%, cukup bukti untuk menyatakan terdapat perbedaan signifikan secara multivariat (Speed dan Accuracy) antara pelatihan X dan pelatihan Y. Terlihat bahwa pelatihan X memberikan nilai akurasi yang lebih tinggi, sedangkan pelatihan Y memberikan nilai speed / kecepatan yang lebih tinggi.

2.4.3 Selang Kepercayaan

2.4.3.1 Simultan

T.ci1 = function(mu_X, mu_Y, S1_1, S2_1, n1, n2, avec=rep(1,length(mu)), level=0.95){
p = length(mu_X)
mu = mu_X-mu_Y
cval = qchisq(level, p)
zhat = crossprod(avec, mu)
zvar = crossprod(avec, S1_1 %*% avec)/n1 + crossprod(avec, S2_1 %*% avec)/n2
const = sqrt(cval * zvar)
c(lower = zhat - const, upper = zhat + const)
}

T.ci1(xbar1, ybar1, S1_1, S2_1, n1, n2, avec=c(1,0),level=0.95) #Speed
##     lower     upper 
## -8.343913 -2.148787
T.ci1(xbar1, ybar1, S1_1, S2_1, n1, n2, avec=c(0,1),level=0.95) #Accuracy
##    lower    upper 
## 2.839469 8.127230

Pada taraf nyata 5%, selang kepercayaan simultan untuk speed berada dalam rentang -8.343 sampai -2.148. Sedangkan selang kepercayaan simultan untuk accuracy berada dalam rentang 2.839 sampai 8.127. Lewat selang kepercayaan simultan ini, dapat disimpulkan bahwa untuk variabel speed maupun accuracy, terdapat perbedaan signifikan antara pelatihan X dengan pelatihan Y.