# ==== S E T U P & D A T A G E N E R A T I O N (NIM 020, Tema: HOTEL) ====
suppressPackageStartupMessages({
library(MASS) # mvrnorm
library(MVTests) # OneSampleHT2
library(Hotelling) # hotelling.test
})
# NIM 020 -> parameter turunan
nim_digits <- c(0,2,0)
set.seed(sum(nim_digits) + 123)
# Ukuran sampel
n_pair <- 24 # berpasangan: 24 tamu menilai dua hotel
n1 <- 20; n2 <- 20 # independen: samakan biar 1 data.frame (A & B) mulus
# Mean & kovarians (Skor Layanan Kamar = Service, Kebersihan = Clean)
# Berpasangan: Hotel A vs Hotel B (tamu yang sama)
mu_A_pair <- c(81, 77) # A_Service, A_Clean
mu_B_pair <- c(76, 74) # B_Service, B_Clean
Sigma_pair <- matrix(c(88, 16,
16, 82), 2, 2)
# Independen (ragam sama)
mu_A_ind <- c(79, 76)
mu_B_ind <- c(83, 80)
Sigma_pool <- matrix(c(132, 18,
18, 126), 2, 2)
# Independen (ragam tidak sama)
Sigma_A <- Sigma_pool
Sigma_B <- matrix(c(148, 21,
21, 140), 2, 2)
# ==== 2. PERBANDINGAN DUA VEKTOR NILAI TENGAH - S A M P E L B E R P A S A N G A N ====
# (Gantikan read_excel sheet=1)
HotelA_pair <- mvrnorm(n_pair, mu = mu_A_pair, Sigma = Sigma_pair)
HotelB_pair <- mvrnorm(n_pair, mu = mu_B_pair, Sigma = Sigma_pair)
# Struktur kolom mengikuti kodenmu: kolom 1-2 kondisi1, kolom 3-4 kondisi2
dataset <- data.frame(
HotelA_Service = HotelA_pair[,1],
HotelA_Clean = HotelA_pair[,2],
HotelB_Service = HotelB_pair[,1],
HotelB_Clean = HotelB_pair[,2]
)
cat("==== [PAIRED] 5 baris pertama dataset ====\n"); print(head(dataset, 5))
## ==== [PAIRED] 5 baris pertama dataset ====
## HotelA_Service HotelA_Clean HotelB_Service HotelB_Clean
## 1 71.00241 74.33882 58.44576 71.14336
## 2 86.98469 78.06244 77.95286 82.09412
## 3 57.20790 77.07546 63.95431 77.03365
## 4 83.33042 72.88337 68.10157 68.65664
## 5 82.65278 68.77271 81.77831 61.34267
# Menghitung selisih (A - B) sesuai indeks
d_kata = dataset[,1] - dataset[,3] # Service
d_kata_kerja = dataset[,2] - dataset[,4] # Clean
X = data.frame(d_kata, d_kata_kerja)
# Rataan & kovarians selisih
xbar = apply(X, 2, mean); cat("\n==== xbar (paired) ====\n"); print(xbar)
##
## ==== xbar (paired) ====
## d_kata d_kata_kerja
## 1.697977 2.806763
cov_m = cov(X); cat("\n==== cov_m (paired) ====\n"); print(cov_m)
##
## ==== cov_m (paired) ====
## d_kata d_kata_kerja
## d_kata 193.11200 70.30478
## d_kata_kerja 70.30478 117.21656
# Uji T^2 Hotelling (One-sample) H0: mu=(0,0)
mean0 = c(0,0)
result = OneSampleHT2(X, mu0 = mean0, alpha = 0.05)
cat("\n==== OneSampleHT2 (paired) ====\n"); print(summary(result))
##
## ==== OneSampleHT2 (paired) ====
## One Sample Hotelling T Square Test
##
## Hotelling T Sqaure Statistic = 1.613031
## F value = 0.771 , df1 = 2 , df2 = 22 , p-value: 0.474
##
## Descriptive Statistics
##
## d_kata d_kata_kerja
## N 24.000000 24.000000
## Means 1.697977 2.806763
## Sd 13.896474 10.826660
##
##
## Detection important variable(s)
##
## Lower Upper Mu0 Important Variables?
## d_kata -5.913301 9.309256 0 FALSE
## d_kata_kerja -3.123138 8.736664 0 FALSE
## NULL
# CI Simultan & Bonferroni (sesuai fungsi kamu)
cat("\n==== CI Simultan (paired) ====\n"); print(result$CI)
##
## ==== CI Simultan (paired) ====
## Lower Upper Mu0 Important Variables?
## d_kata -5.913301 9.309256 0 FALSE
## d_kata_kerja -3.123138 8.736664 0 FALSE
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(X)
cat("\n==== Bonferroni CI (paired) - Service ====\n"); print(bon(xbar, cov_m, n, 0.05, 1))
##
## ==== Bonferroni CI (paired) - Service ====
## lower.d_kata upper.d_kata
## -5.103849 8.499804
cat("\n==== Bonferroni CI (paired) - Clean ====\n"); print(bon(xbar, cov_m, n, 0.05, 2))
##
## ==== Bonferroni CI (paired) - Clean ====
## lower.d_kata_kerja upper.d_kata_kerja
## -2.492500 8.106025
# ==== 3. PERBANDINGAN DUA VEKTOR NILAI TENGAH - I N D E P E N D E N (Ragam Sama) ====
set.seed(sum(nim_digits) + 321) # seed beda agar dataset-nya beda dari paired
HotelA_ind <- mvrnorm(n1, mu = mu_A_ind, Sigma = Sigma_pool)
HotelB_ind <- mvrnorm(n2, mu = mu_B_ind, Sigma = Sigma_pool)
dataset <- data.frame(
HotelA_Service = HotelA_ind[,1],
HotelA_Clean = HotelA_ind[,2],
HotelB_Service = HotelB_ind[,1],
HotelB_Clean = HotelB_ind[,2]
)
cat("\n==== [INDEPENDEN Ragam Sama] 5 baris pertama dataset ====\n"); print(head(dataset, 5))
##
## ==== [INDEPENDEN Ragam Sama] 5 baris pertama dataset ====
## HotelA_Service HotelA_Clean HotelB_Service HotelB_Clean
## 1 75.51586 90.12410 102.16019 89.50635
## 2 68.58600 65.48596 85.76632 68.27814
## 3 71.18833 71.08715 76.69256 81.25085
## 4 77.63853 73.55494 95.30154 73.70658
## 5 69.28173 70.78902 86.47383 74.73001
# Pisah sesuai langkahmu
data_cat = dataset[,1:2] # Hotel A
data_man = dataset[,3:4] # Hotel B
# Rataan, kovarians, & S gabungan
xbar1 = apply(data_cat, 2, mean); cat("\n==== xbar1 (Hotel A) ====\n"); print(xbar1)
##
## ==== xbar1 (Hotel A) ====
## HotelA_Service HotelA_Clean
## 76.95049 74.52897
xbar2 = apply(data_man, 2, mean); cat("\n==== xbar2 (Hotel B) ====\n"); print(xbar2)
##
## ==== xbar2 (Hotel B) ====
## HotelB_Service HotelB_Clean
## 84.94084 77.76504
cov_m1 = cov(data_cat); cat("\n==== cov_m1 (A) ====\n"); print(cov_m1)
##
## ==== cov_m1 (A) ====
## HotelA_Service HotelA_Clean
## HotelA_Service 72.27193 -36.24243
## HotelA_Clean -36.24243 151.81119
cov_m2 = cov(data_man); cat("\n==== cov_m2 (B) ====\n"); print(cov_m2)
##
## ==== cov_m2 (B) ====
## HotelB_Service HotelB_Clean
## HotelB_Service 199.59626 62.92135
## HotelB_Clean 62.92135 136.47373
n1 = nrow(data_cat); n2 = nrow(data_man)
s_gab = ((n1-1)*cov_m1 + (n2-1)*cov_m2) / (n1 + n2 - 2)
cat("\n==== S Gabungan ====\n"); print(s_gab)
##
## ==== S Gabungan ====
## HotelA_Service HotelA_Clean
## HotelA_Service 135.93409 13.33946
## HotelA_Clean 13.33946 144.14246
# Hotelling T^2 dua sampel (var.equal = TRUE)
t2_homogen = hotelling.test(data_cat, data_man, var.equal = TRUE)
cat("\n==== Hotelling T^2 (var.equal=TRUE) ====\n"); print(t2_homogen)
##
## ==== Hotelling T^2 (var.equal=TRUE) ====
## Test stat: 5.1177
## Numerator df: 2
## Denominator df: 37
## P-value: 0.09657
# CI Simultan & Bonferroni
T.ci = function(mu1, mu2, S_gab, n1, n2, avec = rep(1, length(mu1)), level = 0.95){
p = length(mu1)
mu = mu1 - mu2
cval = qf(level, p, n1 + n2 - p - 1) * p * (n1 + n2 - 2) / (n1 + n2 - p - 1)
zhat = as.numeric(crossprod(avec, mu))
zvar = as.numeric(crossprod(avec, S_gab %*% avec)) * (1/n1 + 1/n2)
const = sqrt(cval * zvar)
c(lower = zhat - const, upper = zhat + const)
}
cat("\n==== CI Simultan (Ragam Sama) - Service ====\n"); print(T.ci(xbar1, xbar2, s_gab, n1, n2, avec = c(1,0), level = 0.95))
##
## ==== CI Simultan (Ragam Sama) - Service ====
## lower upper
## -17.519198 1.538494
cat("\n==== CI Simultan (Ragam Sama) - Clean ====\n"); print(T.ci(xbar1, xbar2, s_gab, n1, n2, avec = c(0,1), level = 0.95))
##
## ==== CI Simultan (Ragam Sama) - Clean ====
## lower upper
## -13.048391 6.576266
bon_ind = function(mu1, mu2, S, n1, n2, alpha, k){
p = length(mu1)
mu = mu1 - mu2
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
}
cat("\n==== Bonferroni CI (Ragam Sama) - Service ====\n"); print(bon_ind(xbar1, xbar2, s_gab, n1, n2, 0.05, 1))
##
## ==== Bonferroni CI (Ragam Sama) - Service ====
## lower upper
## -16.5946043 0.6139007
cat("\n==== Bonferroni CI (Ragam Sama) - Clean ====\n"); print(bon_ind(xbar1, xbar2, s_gab, n1, n2, 0.05, 2))
##
## ==== Bonferroni CI (Ragam Sama) - Clean ====
## lower upper
## -12.096291 5.624166
# ==== 4. PERBANDINGAN DUA VEKTOR NILAI TENGAH - I N D E P E N D E N (Ragam Tidak Sama) ====
set.seed(sum(nim_digits) + 456) # seed berbeda agar dataset berbeda dari ragam sama
HotelA_uneq <- mvrnorm(n1, mu = mu_A_ind, Sigma = Sigma_A)
HotelB_uneq <- mvrnorm(n2, mu = mu_B_ind, Sigma = Sigma_B)
dataset <- data.frame(
HotelA_Service = HotelA_uneq[,1],
HotelA_Clean = HotelA_uneq[,2],
HotelB_Service = HotelB_uneq[,1],
HotelB_Clean = HotelB_uneq[,2]
)
cat("\n==== [INDEPENDEN Ragam Tidak Sama] 5 baris pertama dataset ====\n"); print(head(dataset, 5))
##
## ==== [INDEPENDEN Ragam Tidak Sama] 5 baris pertama dataset ====
## HotelA_Service HotelA_Clean HotelB_Service HotelB_Clean
## 1 84.12128 84.12258 94.22681 73.94576
## 2 84.63697 74.07498 85.05468 80.22862
## 3 73.26890 64.48100 80.86816 77.80220
## 4 78.75404 79.74370 115.95055 97.04098
## 5 55.84222 81.07095 78.34920 90.53228
data_cat = dataset[,1:2]
data_man = dataset[,3:4]
xbar1 = apply(data_cat, 2, mean); cov_m1 = cov(data_cat)
xbar2 = apply(data_man, 2, mean); cov_m2 = cov(data_man)
n1 = nrow(data_cat); n2 = nrow(data_man)
# Hotelling T^2 dua sampel (var.equal = FALSE)
t2_not_homogen = hotelling.test(data_cat, data_man, var.equal = FALSE)
cat("\n==== Hotelling T^2 (var.equal=FALSE) ====\n"); print(t2_not_homogen)
##
## ==== Hotelling T^2 (var.equal=FALSE) ====
## Test stat: 1.9817
## Numerator df: 2
## Denominator df: 35.8056657814176
## P-value: 0.3908
# CI Simultan (chi-square)
T.ci.uneq = function(mu1, mu2, S1, S2, n1, n2, avec = rep(1, length(mu1)), level = 0.95){
p = length(mu1)
mu = mu1 - mu2
cval = qchisq(level, p)
zhat = as.numeric(crossprod(avec, mu))
zvar = as.numeric(crossprod(avec, S1 %*% avec))/n1 + as.numeric(crossprod(avec, S2 %*% avec))/n2
const = sqrt(cval * zvar)
c(lower = zhat - const, upper = zhat + const)
}
cat("\n==== CI Simultan (Ragam Tidak Sama) - Service ====\n"); print(T.ci.uneq(xbar1, xbar2, cov_m1, cov_m2, n1, n2, avec = c(1,0), level = 0.95))
##
## ==== CI Simultan (Ragam Tidak Sama) - Service ====
## lower upper
## -14.481550 5.127496
cat("\n==== CI Simultan (Ragam Tidak Sama) - Clean ====\n"); print(T.ci.uneq(xbar1, xbar2, cov_m1, cov_m2, n1, n2, avec = c(0,1), level = 0.95))
##
## ==== CI Simultan (Ragam Tidak Sama) - Clean ====
## lower upper
## -9.927601 4.201285