set.seed(012)
library(MASS)
## Warning: package 'MASS' was built under R version 4.4.3
n1 <- 20
n2 <- 20
mu_A <- c(5, 30)
mu_B <- c(7, 28)
Sigma <- matrix(c(16, 4,
4, 20), 2, 2, byrow = TRUE)
A <- mvrnorm(n1, mu = mu_A, Sigma = Sigma)
B <- mvrnorm(n2, mu = mu_B, Sigma = Sigma)
colnames(A) <- colnames(B) <- c("Sales","Profit")
head(A); head(B)
## Sales Profit
## [1,] 0.6103941 24.46207
## [2,] 2.6506883 40.24115
## [3,] -0.5506086 28.09875
## [4,] 3.6534559 25.70524
## [5,] 3.2291408 19.96206
## [6,] 5.1579500 28.38494
## Sales Profit
## [1,] 5.989675 29.20947
## [2,] 1.007757 25.26133
## [3,] 5.763382 31.98608
## [4,] 2.406595 21.94683
## [5,] 3.620777 28.36926
## [6,] 5.470802 31.44987
xbar1 <- colSums(A) / n1
xbar2 <- colSums(B) / n2
mu_hat <- xbar1 - xbar2
xbar1; xbar2; mu_hat
## Sales Profit
## 3.868968 28.853266
## Sales Profit
## 6.365435 28.063104
## Sales Profit
## -2.496467 0.790162
center1 <- scale(A, center = xbar1, scale = FALSE)
S1 <- t(center1) %*% center1 / (n1 - 1)
center2 <- scale(B, center = xbar2, scale = FALSE)
S2 <- t(center2) %*% center2 / (n2 - 1)
S1; S2
## Sales Profit
## Sales 9.446645 1.389631
## Profit 1.389631 18.011733
## Sales Profit
## Sales 16.58404 3.98957
## Profit 3.98957 12.19236
Sp <- ((n1-1)*S1 + (n2-1)*S2) / (n1 + n2 - 2)
Sp
## Sales Profit
## Sales 13.01534 2.68960
## Profit 2.68960 15.10205
p <- 2
T2 <- (n1*n2/(n1+n2)) * t(mu_hat) %*% solve(Sp) %*% mu_hat
Fh <- as.numeric(((n1+n2-p-1)/(p*(n1+n2-2))) * T2)
T2; Fh
## [,1]
## [1,] 5.961119
## [1] 2.902124
Fcrit <- qf(0.95, df1 = 2, df2 = (n1+n2-2))
pval <- 1 - pf(Fh, df1 = 2, df2 = (n1+n2-2))
Fcrit; pval
## [1] 3.244818
## [1] 0.06715526
Hasil uji menunjukkan bahwa nilai Fhitung lebih kecil daripada \(F_{kritis}\) =3.244818 , dan nilai \(p-value\) yang diperoleh 0.0672 >0.05
Dengan demikian, tak tolak H0 pada taraf signifikansi 5%.
Artinya, secara multivariat (Sales dan Profit) belum terdapat cukup bukti untuk mengatakan antara grup A (model lama) dan grup B (model baru) berbeda signifikan.
set.seed(012)
library(MASS)
n1 <- 20
n2 <- 20
mu_A <- c(5, 30)
mu_B <- c(7, 28)
Sigma <- matrix(c(16, 4,
4, 20), 2, 2, byrow = TRUE)
A <- mvrnorm(n1, mu = mu_A, Sigma = Sigma)
B <- mvrnorm(n2, mu = mu_B, Sigma = Sigma)
colnames(A) <- colnames(B) <- c("Sales","Profit")
data_sales = A # Toko A
data_profit = B # Toko B
xbar1 = apply(data_sales, 2, mean)
xbar1
## Sales Profit
## 3.868968 28.853266
xbar2 = apply(data_profit, 2, mean)
xbar2
## Sales Profit
## 6.365435 28.063104
cov_m1 = cov(data_sales)
cov_m1
## Sales Profit
## Sales 9.446645 1.389631
## Profit 1.389631 18.011733
cov_m2 = cov(data_profit)
cov_m2
## Sales Profit
## Sales 16.58404 3.98957
## Profit 3.98957 12.19236
n1 = nrow(data_sales)
n2 = nrow(data_profit)
s_gab = ((n1-1)*cov_m1+(n2-1)*cov_m2)/(n1+n2-2)
s_gab
## Sales Profit
## Sales 13.01534 2.68960
## Profit 2.68960 15.10205
library(Hotelling)
## Loading required package: corpcor
t2_homogen = hotelling.test(data_sales, data_profit, var.equal=TRUE)
t2_homogen
## Test stat: 5.9611
## Numerator df: 2
## Denominator df: 37
## P-value: 0.06749
T.ci = function(mu1, mu2, S_gab, n1, n2, avec=rep(1,length(mu)), 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 = crossprod(avec, mu)
zvar = crossprod(avec, S_gab %*% avec)* (1/n1+1/n2)
const = sqrt(cval * zvar)
c(lower = zhat - const, upper = zhat + const)
}
# Sales
T.ci(xbar1, xbar2, s_gab, n1, n2, avec=c(1,0), level=0.95)
## lower upper
## -5.4449859 0.4520517
# Profit
T.ci(xbar1, xbar2, s_gab, n1, n2, avec=c(0,1), level=0.95)
## lower upper
## -2.385937 3.966261
bon = 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
}
# Sales
bon(xbar1, xbar2, s_gab, n1, n2, 0.05, 1)
## lower upper
## -5.1588882 0.1659539
# Profit
bon(xbar1, xbar2, s_gab, n1, n2, 0.05, 2)
## lower upper
## -2.077757 3.658081
Hasil uji menunjukkan bahwa nilai p-value yang diperoleh \(p-value\) 0.06749 >0.05. Dengan demikian, tak tolak H0 pada taraf signifikansi 5%. Artinya, Sales dan Profit antara Toko A dan Toko B belum terdapat cukup bukti untuk mengatakan antara Toko A dan Toko B berbeda signifikan.
library(MASS)
library(Hotelling)
set.seed(012)
n_online <- 50
n_offline <- 50
# --- sementara kita simulasikan data awal agar contoh bisa jalan ---
# (kalau sudah punya data asli, lewati bagian ini)
mu_online_true <- c(65, 70)
mu_offline_true <- c(75, 85)
Sigma_true <- matrix(c(150, 30,
30, 170), nrow = 2)
data_online <- mvrnorm(n_online, mu = mu_online_true, Sigma = Sigma_true)
data_offline <- mvrnorm(n_offline, mu = mu_offline_true, Sigma = Sigma_true)
data_mahasiswa <- data.frame(
Platform = rep(c("Online","Offline"), each = 50),
Pengetahuan = c(data_online[,1], data_offline[,1]),
Analitik = c(data_online[,2], data_offline[,2])
)
online_data <- subset(data_mahasiswa, Platform == "Online",
select = c(Pengetahuan, Analitik))
offline_data <- subset(data_mahasiswa, Platform == "Offline",
select = c(Pengetahuan, Analitik))
# Rata-rata (μ) otomatis
mu_online <- colMeans(online_data)
mu_offline <- colMeans(offline_data)
# Matriks kovarians (Σ) otomatis
Sigma_online <- cov(online_data)
Sigma_offline <- cov(offline_data)
# Jika ingin kovarians gabungan (pooled)
Sigma_pooled <- cov(rbind(online_data, offline_data))
# Cek hasil
mu_online
## Pengetahuan Analitik
## 63.10235 68.92884
mu_offline
## Pengetahuan Analitik
## 75.50575 85.31064
Sigma_online
## Pengetahuan Analitik
## Pengetahuan 121.52686 26.80425
## Analitik 26.80425 116.79474
Sigma_offline
## Pengetahuan Analitik
## Pengetahuan 138.04757 40.08765
## Analitik 40.08765 190.03140
Sigma_pooled
## Pengetahuan Analitik
## Pengetahuan 167.32584 84.41872
## Analitik 84.41872 219.63191
dataset <- data.frame(
Online_Pengetahuan = online_data$Pengetahuan,
Online_Analitik = online_data$Analitik,
Offline_Pengetahuan = offline_data$Pengetahuan,
Offline_Analitik = offline_data$Analitik
)
t2_not_homogen <- hotelling.test(
dataset[,1:2], # Online
dataset[,3:4], # Offline
var.equal = FALSE
)
t2_not_homogen
## Test stat: 59.652
## Numerator df: 2
## Denominator df: 95.2760326274497
## P-value: 1.054e-10
n1 <- n_online
n2 <- n_offline
S_gab <- ((n1 - 1) * Sigma_online + (n2 - 1) * Sigma_offline) / (n1 + n2 - 2)
# Untuk variabel Pengetahuan
T.ci(mu_online, mu_offline, S_gab, n1, n2,
avec = c(1, 0), level = 0.95)
## lower upper
## -18.096919 -6.709889
# Untuk variabel Analitik
T.ci(mu_online, mu_offline, S_gab, n1, n2,
avec = c(0, 1), level = 0.95)
## lower upper
## -22.57186 -10.19172
Dengan nilai \(p-value\) yang sangat kecil (<0.05), maka tolak \(H_0\). Terdapat cukup bukti untuk mengatakan bahwa terdapat perbedaan antara nilai mahasiswa yang mengikuti pelajaran secara online berbeda signifikan dengan nilai mahasiswa yang mengikuti pelajaran secara offline.