Sebuah perusahaan ingin mengevaluasi efektivitas program pelatihan produktivitas karyawan. Dari 25 karyawan, diukur dua hal:
Produktivitas kerja (jumlah unit yang diselesaikan per minggu),
Kepuasan kerja (skor skala 0–100).
Pengukuran dilakukan sebelum pelatihan dan setelah pelatihan pada orang yang sama. Karena datanya berasal dari orang yang sama dalam dua kondisi berbeda, maka sampel berpasangan.
set.seed(028)
library(MASS)## Warning: package 'MASS' was built under R version 4.4.3n <- 25
mu_before <- c(70, 65)
mu_after <- c(78, 72)
Sigma <- matrix(c(20,10,10,30), nrow=2)
# Simulasi data
before <- mvrnorm(n, mu_before, Sigma)
after <- mvrnorm(n, mu_after, Sigma)
df1 <- data.frame(
Prod_before = before[,1],
Kepuasan_before = before[,2],
Prod_after = after[,1],
Kepuasan_after = after[,2])
head(df1)## Prod_before Kepuasan_before Prod_after Kepuasan_after
## 1 65.23110 54.49706 83.32936 73.37639
## 2 67.82721 65.88822 69.66162 71.45060
## 3 71.04051 54.94416 71.16619 78.70148
## 4 65.14686 55.13012 73.63777 64.65114
## 5 73.30652 64.10671 70.16909 60.17894
## 6 74.72497 65.83735 76.65106 71.36631xbar = apply(df1, 2, mean)
xbar## Prod_before Kepuasan_before Prod_after Kepuasan_after
## 69.62211 64.11534 75.75710 71.22803cov_m = cov(df1)
cov_m## Prod_before Kepuasan_before Prod_after Kepuasan_after
## Prod_before 20.851819 9.941594 1.607731 2.035949
## Kepuasan_before 9.941594 33.411105 1.230634 2.005932
## Prod_after 1.607731 1.230634 15.112399 5.747214
## Kepuasan_after 2.035949 2.005932 5.747214 26.248556library(MVTests)## Warning: package 'MVTests' was built under R version 4.4.3##
## Attaching package: 'MVTests'## The following object is masked from 'package:datasets':
##
## irismean0 <- colMeans(df1) # rata-rata tiap kolom Kepuasan dan Produksi
result <- OneSampleHT2(df1, mu0 = mean0, alpha = 0.05)
summary(result)## One Sample Hotelling T Square Test
##
## Hotelling T Sqaure Statistic = 0
## F value = 0 , df1 = 4 , df2 = 21 , p-value: 1
##
## Descriptive Statistics
##
## Prod_before Kepuasan_before Prod_after Kepuasan_after
## N 25.000000 25.000000 25.000000 25.000000
## Means 69.622113 64.115344 75.757097 71.228033
## Sd 4.566379 5.780234 3.887467 5.123334
##
##
## Detection important variable(s)
##
## Lower Upper Mu0 Important Variables?
## Prod_before 66.33136 72.91286 69.62211 FALSE
## Kepuasan_before 59.94983 68.28085 64.11534 FALSE
## Prod_after 72.95560 78.55859 75.75710 FALSE
## Kepuasan_after 67.53592 74.92015 71.22803 FALSEKarena statistik uji = 0 dan p-value = 1, maka tidak ada bukti sama sekali untuk menolak hipotesis nol (H₀). Artinya, vektor rata-rata selisih (antara kondisi sesudah dan sebelum pelatihan) hampir persis nol dalam data yang diuji. Dengan kata lain, secara multivariat (mencakup semua variabel yang dimasukkan sekaligus), tidak ada perbedaan yang signifikan antara sebelum dan sesudah.
result$CI## Lower Upper Mu0 Important Variables?
## Prod_before 66.33136 72.91286 69.62211 FALSE
## Kepuasan_before 59.94983 68.28085 64.11534 FALSE
## Prod_after 72.95560 78.55859 75.75710 FALSE
## Kepuasan_after 67.53592 74.92015 71.22803 FALSEHasil selang kepercayaan simultan menunjukkan bahwa baik produksi maupun kepuasan karyawan mengalami peningkatan setelah pelatihan. Rentang rata-rata setelah pelatihan bergeser ke atas dibandingkan sebelum pelatihan, sehingga dapat disimpulkan bahwa pelatihan memberi dampak positif terhadap kinerja dan kepuasan kerja karyawan.
bon = function(mu,S,n,alpha,k){
p = length(mu)
lower = mu[k] - sqrt(S[k,k]/n) * abs(qt(alpha/(2*p), df=n-1))
upper = mu[k] + sqrt(S[k,k]/n) * abs(qt(alpha/(2*p), df=n-1))
c(lower = lower,upper = upper)
}
n = nrow(df)#Kepuasan
bon(xbar, cov_m,n,0.05,3)## numeric(0)Hasil selang kepercayaan Bonferroni untuk variabel produktivitas adalah dari 73.65 hingga 77.856.Karena interval ini tidak mencakup nol, maka dapat disimpulkan bahwa terdapat perbedaan yang signifikan pada produktivitas karyawan sebelum dan sesudah pelatihan. Artinya, program pelatihan memang memberikan peningkatan produktivitas yang nyata pada taraf kepercayaan 95%.
#Produksi
bon(xbar, cov_m,n,0.05,4)## numeric(0)Hasil selang kepercayaan Bonferroni untuk kepuasan sesudah pelatihan adalah 68.46 – 73.99. Ini menunjukkan bahwa rata-rata kepuasan kerja karyawan cenderung lebih tinggi setelah pelatihan.Dengan kata lain, selain meningkatkan produktivitas, pelatihan juga memberi dampak positif terhadap kepuasan karyawan.
Perbedaan hasil uji T2 Hotelling dengan hasil Selang Kepercayaan Simultan maupun Bonferroni adalah Hotelling melihat perbedaan rata-rata secara keseluruhan (multivariat), sedangkan SK simultan memeriksa satu per satu variabel dengan koreksi. Karena itu hasilnya bisa berbeda.
Sebuah perusahaan ingin membandingkan dua populasi karyawan dari dua cabang berbeda mengenai kinerja kerja mereka.
Variabel yang diukur ada dua (multivariat):
Produktivitas kerja (jumlah unit yang dihasilkan per minggu).
Tingkat kepuasan kerja (skor survei skala 1–100).
# ===== Persiapan =====
set.seed(028) # agar reproduktif
library(MASS) # untuk mvrnorm
# ===== Parameter populasi =====
# Means (nilai tengah masing-masing variabel)
mu_A <- c(Produktivitas = 120, Kepuasan = 75) # Cabang A
mu_B <- c(Produktivitas = 110, Kepuasan = 68) # Cabang B (sedikit lebih rendah)
# Matriks kovarian sama untuk kedua populasi (2x2)
# Var(Produktivitas) = 100, Var(Kepuasan) = 225, Cov = 50
Sigma <- matrix(c(100, 50,
50, 225), nrow = 2, byrow = TRUE)
# Ukuran sampel
nA <- 30
nB <- 35
# ===== Membangkitkan data masing-masing grup =====
data_A <- mvrnorm(n = nA, mu = mu_A, Sigma = Sigma)
data_B <- mvrnorm(n = nB, mu = mu_B, Sigma = Sigma)
# Ubah jadi data.frame dan tambahkan label grup
df_A <- data.frame(data_A)
df_A$Group <- "A"
df_B <- data.frame(data_B)
df_B$Group <- "B"
# Gabungkan
df <- rbind(df_A, df_B)
rownames(df) <- NULL
# ===== Ambil ukuran sampel =====
n1 <- nrow(df_A)
n2 <- nrow(df_B)
# ===== Hitung kovarian masing-masing =====
cov_m1 <- cov(df_A[, 1:2]) # hanya ambil variabel Produktivitas & Kepuasan
cov_m2 <- cov(df_B[, 1:2])
# ===== Rumus pooled covariance matrix =====
s_gab <- ((n1 - 1) * cov_m1 + (n2 - 1) * cov_m2) / (n1 + n2 - 2)
s_gab## Produktivitas Kepuasan
## Produktivitas 114.86590 38.06044
## Kepuasan 38.06044 222.65556Uji T2 Hotelling Dua Populasi Sampel Saling Bebas Ragam Sama
library(Hotelling)## Loading required package: corpcort2_homogen = hotelling.test(data_A,data_B,var.equal=TRUE)
t2_homogen## Test stat: 4.821
## Numerator df: 2
## Denominator df: 62
## P-value: 0.1017Berdasarkan hasil uji Hotelling’s T², diperoleh nilai statistik uji sebesar 4.821 dengan derajat bebas numerator = 2 serta p-value = 0.1017. Karena p-value lebih besar dari taraf signifikansi 5%, maka keputusan uji adalah gagal menolak H0.
T.ci = function(mu_A, mu_B, S_gab, n1, n2, avec=rep(1,length(mu)), level=0.95){
p = length(mu_A)
mu = mu_A-mu_B
cval = qf(level, p, n1+n2-p-1) * p * (n1+n2-2) / (n1+n2-p-1)
zhat = crossprod(avec, mu)
zvar = crossprod(avec, S_gab %*% avec)* (1/n1+1/n2)
const = sqrt(cval * zvar)
c(lower = zhat - const, upper = zhat + const)
}#Produktivitas
xbar1 = apply(data_A, 2, mean)
xbar1## Produktivitas Kepuasan
## 119.34406 71.62496xbar2 = apply(data_B, 2, mean)
xbar2## Produktivitas Kepuasan
## 113.51133 69.00266T.ci(xbar1, xbar2, s_gab, n1,n2, avec=c(1,0),level=0.95)## lower upper
## -0.9090609 12.5745176Hasil uji simultanmenghasilkan nilai statistik sebesar 4.821 dengan p-value 0.1017. Karena p-value lebih besar dari taraf signifikansi 5%, maka keputusan uji adalah gagal menolak H0. Artinya, secara simultan tidak terdapat perbedaan yang signifikan pada rata-rata multivariat (Produktivitas dan Kepuasan) antara Cabang A dan Cabang B. Dengan kata lain, kedua cabang dapat dianggap memiliki profil produktivitas dan kepuasan yang relatif sama jika dilihat bersama-sama.
#Kepuasan
T.ci(xbar1, xbar2, s_gab, n1,n2, avec=c(0,1),level=0.95)## lower upper
## -6.764041 12.008643Interval kepercayaan simultan Bonferroni Kepuasan berada pada rentang 6.76 hingga 12.01. Karena rentang ini masih mencakup nilai nol, maka secara statistik tidak terdapat bukti yang cukup untuk menyatakan adanya perbedaan signifikan tingkat kepuasan antara Cabang A dan Cabang B. Dengan demikian, meskipun rata-rata kepuasan Cabang A terlihat sedikit lebih tinggi, secara simultan perbedaan tersebut tidak signifikan.
bon = function(mu_A, mu_B ,S, n1, n2, alpha, k){
p = length(mu_A)
mu = mu_A-mu_B
lower = mu[k] - sqrt((S[k,k]) *(1/n1+1/n2))* abs(qt(alpha/(2*p), df=n1+n2-2))
upper = mu[k] + sqrt((S[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
}#Produktivitas
bon(xbar1, xbar2, s_gab, n1, n2,0.05,1)## lower upper
## -0.2904125 11.9558692Hasil interval kepercayaan Bonferroni untuk variabel Produktivitas berada pada rentang 0.29 hingga 11.96.Namun karena interval tersebut masih mencakup angka nol, maka secara statistik pada taraf signifikansi yang sudah disesuaikan dengan Bonferroni, perbedaan rata-rata Produktivitas belum signifikan.
#Kepuasan
bon(xbar1, xbar2, s_gab, n1, n2,0.05,2)## lower upper
## -5.90272 11.14732Hasil interval kepercayaan Bonferroni diperoleh dari 5.90 hingga 11.15. Selisih rata-rata antara Cabang A dan Cabang B berada di sekitar 2.62 poin, tetapi karena rentang interval masih mencakup nol, maka secara statistik tidak ada perbedaan signifikan antara tingkat kepuasan karyawan di kedua cabang.
Sebuah rumah sakit membandingkan dua unit rehabilitasi pasca-operasi (Unit X vs Unit Y) pada dua indikator pasien: Waktu Pemulihan (hari) dan Skor Fungsi Harian (0–100). Unit X menerapkan program pemulihan individual yang hasilnya sangat bervariasi antar pasien (beberapa cepat, beberapa lambat) — sehingga varians lebih besar. Unit Y memakai protokol standar ketat sehingga hasilnya lebih konsisten (varians lebih kecil).
set.seed(028)
# Parameter
mu_X <- c(Waktu = 18, Fungsi = 70)
mu_Y <- c(Waktu = 15, Fungsi = 66)
Sigma_X <- matrix(c(40, -8, -8, 64), nrow=2)
Sigma_Y <- matrix(c(15, 2, 2, 30), nrow=2)
nX <- 40; nY <- 30
# Simulasi
data_X <- mvrnorm(n = nX, mu = mu_X, Sigma = Sigma_X)
data_Y <- mvrnorm(n = nY, mu = mu_Y, Sigma = Sigma_Y)
df_X <- data.frame(data_X); df_X$Group <- "X"
df_Y <- data.frame(data_Y); df_Y$Group <- "Y"
df <- rbind(df_X, df_Y)
colnames(df)[1:2] <- c("Waktu","Fungsi")
df## Waktu Fungsi Group
## 1 19.342141 54.20882 X
## 2 16.489148 68.99506 X
## 3 19.348089 59.07280 X
## 4 22.580858 55.88909 X
## 5 15.507092 70.63040 X
## 6 24.870323 76.60520 X
## 7 4.473121 68.42084 X
## 8 30.269290 73.89037 X
## 9 24.113084 69.20035 X
## 10 11.649299 83.75804 X
## 11 15.674731 63.70864 X
## 12 19.141055 83.29195 X
## 13 15.816536 69.79839 X
## 14 28.122494 66.61805 X
## 15 36.213978 59.68612 X
## 16 16.668762 78.78664 X
## 17 18.905613 81.81484 X
## 18 18.733606 66.51415 X
## 19 6.857302 67.84386 X
## 20 21.151355 74.22431 X
## 21 20.497756 69.11379 X
## 22 16.416601 64.22998 X
## 23 23.520903 66.94719 X
## 24 19.269211 76.49896 X
## 25 32.299749 61.09986 X
## 26 21.503585 67.70490 X
## 27 22.370075 76.62653 X
## 28 24.025043 57.68699 X
## 29 16.737401 67.18835 X
## 30 19.578937 62.95948 X
## 31 17.912065 61.77503 X
## 32 18.275632 85.23315 X
## 33 14.613128 72.05197 X
## 34 10.329254 81.55298 X
## 35 23.095696 70.81035 X
## 36 21.237470 80.98671 X
## 37 9.039101 60.87719 X
## 38 5.202965 59.13816 X
## 39 17.384290 72.90300 X
## 40 18.750964 64.97801 X
## 41 23.889665 66.05040 Y
## 42 15.831374 61.35092 Y
## 43 17.973354 71.88126 Y
## 44 22.871918 74.84740 Y
## 45 23.195542 76.09051 Y
## 46 12.302159 70.67579 Y
## 47 24.470272 63.12078 Y
## 48 11.144266 64.55221 Y
## 49 16.980889 58.32515 Y
## 50 17.221069 70.02314 Y
## 51 17.021465 66.37860 Y
## 52 18.962395 63.98037 Y
## 53 13.109806 68.97005 Y
## 54 16.549132 70.16474 Y
## 55 10.275368 70.35201 Y
## 56 17.111717 71.15175 Y
## 57 20.913219 67.47989 Y
## 58 6.878468 67.14388 Y
## 59 14.654626 61.05032 Y
## 60 25.382215 67.48866 Y
## 61 13.895855 70.20420 Y
## 62 18.840018 68.89371 Y
## 63 11.830182 68.84419 Y
## 64 18.818556 62.83195 Y
## 65 19.004881 62.66492 Y
## 66 20.772121 65.75935 Y
## 67 18.511962 67.08432 Y
## 68 16.537821 65.23332 Y
## 69 10.226650 55.30309 Y
## 70 23.175260 62.11374 Yxbar1 = apply(data_X, 2, mean)
xbar1## Waktu Fungsi
## 18.94969 69.33301xbar2 = apply(data_Y, 2, mean)
xbar2## Waktu Fungsi
## 17.27841 66.66702cov_m1 = cov(data_X)
cov_m1## Waktu Fungsi
## Waktu 43.776353 -6.068955
## Fungsi -6.068955 66.793236cov_m2 = cov(data_Y)
cov_m2## Waktu Fungsi
## Waktu 22.010041 2.528531
## Fungsi 2.528531 21.483450n1 = nrow(data_X)
n2 = nrow(data_Y)library(Hotelling)
t2_not_homogen = hotelling.test(data_X,data_Y,var.equal=FALSE)
t2_not_homogen## Test stat: 4.6495
## Numerator df: 2
## Denominator df: 65.6062061004143
## P-value: 0.1092Berdasarkan hasil uji Hotelling’s T², diperoleh nilai statistik uji sebesar 4.6495 dengan derajat bebas numerator = 2, serta nilai p-value = 0.1092. Karena p-value lebih besar dari taraf signifikansi, maka keputusan uji adalah gagal menolak hipotesis nol (H₀). Artinya, secara simultan tidak terdapat bukti yang cukup untuk menyimpulkan adanya perbedaan rata-rata multivariat (waktu pemulihan dan skor fungsi harian) antara Unit X dan Unit Y. Perbedaan yang muncul lebih mungkin disebabkan oleh variasi sampel daripada benar-benar mencerminkan perbedaan populasi.
T.ci = function(mu1, mu2, S1, S2, n1, n2, avec=rep(1,length(mu)), level=0.95){
p = length(mu1)
mu = mu1-mu2
cval = qchisq(level, p)
zhat = crossprod(avec, mu)
zvar = crossprod(avec, S1 %*% avec)/n1 + crossprod(avec, S2 %*% avec)/n2
const = sqrt(cval * zvar)
c(lower = zhat - const, upper = zhat + const)
}
#Waktu Pemulihan
T.ci(xbar1, xbar2, cov_m1, cov_m2, n1,n2, avec=c(1,0),level=0.95)## lower upper
## -1.638225 4.980795Karena selang ini mencakup nilai nol, maka tidak ada bukti signifikan bahwa rata-rata waktu pemulihan antara Unit X dan Unit Y berbeda secara nyata pada taraf kepercayaan 95%.
#SKOR Harian
T.ci(xbar1, xbar2, cov_m1, cov_m2, n1,n2, avec=c(0,1),level=0.95)## lower upper
## -1.114923 6.446905