#Menyiapkan Dataset
library(readxl)
dataset_4 <- read_excel("F:/APG_Praktikum/APG_Bu Sarni/Pertemuan 4/dataset_4.xlsx")
head(dataset_4,10)## # A tibble: 10 x 8
## X11 X21 X31 X41 X12 X22 X32 X42
## <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
## 1 15 17 24 14 13 14 12 21
## 2 17 15 32 26 14 12 14 26
## 3 15 14 29 23 12 19 21 21
## 4 13 12 10 16 12 13 10 16
## 5 20 17 26 28 11 20 16 16
## 6 15 21 26 21 12 9 14 18
## 7 15 13 26 22 10 13 18 24
## 8 13 5 22 22 10 8 13 23
## 9 14 7 30 17 12 20 19 23
## 10 17 15 30 27 11 10 11 27populasi1 <- matrix(c(dataset_4$X11,dataset_4$X21,dataset_4$X31,dataset_4$X41),32,4)
populasi2 <- matrix(c(dataset_4$X12,dataset_4$X22,dataset_4$X32,dataset_4$X42),32,4)
#Modifikasi
gabung=rbind(populasi1,populasi2)
colnames(gabung)=c("x1","x2","x3","x4")
ind_populasi=factor(c(rep(1,32),rep(2,32)))#Analisis Deskriptif
summary(dataset_4)## X11 X21 X31 X41
## Min. :10.00 Min. : 5.00 Min. :10.00 Min. :14.00
## 1st Qu.:15.00 1st Qu.:14.00 1st Qu.:25.75 1st Qu.:21.00
## Median :16.00 Median :16.00 Median :28.50 Median :23.00
## Mean :15.97 Mean :15.91 Mean :27.19 Mean :22.75
## 3rd Qu.:18.00 3rd Qu.:18.25 3rd Qu.:30.25 3rd Qu.:26.00
## Max. :20.00 Max. :21.00 Max. :34.00 Max. :29.00
## X12 X22 X32 X42
## Min. : 2.00 Min. : 5.00 Min. : 7.00 Min. : 9.00
## 1st Qu.:11.75 1st Qu.:10.00 1st Qu.:13.00 1st Qu.:18.00
## Median :13.00 Median :15.00 Median :16.00 Median :23.00
## Mean :12.34 Mean :13.91 Mean :16.72 Mean :21.94
## 3rd Qu.:14.00 3rd Qu.:17.00 3rd Qu.:21.25 3rd Qu.:26.00
## Max. :17.00 Max. :21.00 Max. :26.00 Max. :28.00#Asumsi Sigma1 = Sigma2
multivariate<-function(mdt1, mdt2,Ho,alfa=0.05){
#populasi 1
n1<-nrow(mdt1)
p1<-ncol(mdt1)
#Membentuk vektor rata-rata sampel
xbar1=apply(mdt1, 2, mean)
varcovar1=cov(mdt1)
#populasi 2
n2<-nrow(mdt2)
p2<-ncol(mdt2)
#Membentuk vektor rata-rata sampel
xbar2=apply(mdt2, 2, mean)
varcovar2=cov(mdt2)
#membentuk S_pooled
spl=((n1-1)/(n1+n2-2))*varcovar1 + ((n2-1)/(n1+n2-2))*varcovar2
#T hotelling
Tsq<-(n1*n2/(n1+n2))*(t((xbar1-xbar2)-Ho)%*%solve(spl)%*%((xbar1-xbar2)-Ho))
#wilayah kritis
Ftab<-qf(1-alfa,p1,n1+n2-p1-1)
csq<-((n1+n2-2)*p1/(n1+n2-p1-1))*Ftab
uji<-ifelse(Tsq>csq, "Tolak Ho","Gagal Tolak Ho")
hasil<-list(xbar1, varcovar1, xbar2, varcovar2,spl,Tsq,Ftab,csq,uji)
names(hasil)<-c("Vektor rata 1","matriks var-covar1","Vektor rata 2",
"matriks var-covar2","s-pooled","T-square","F tabel","c-square","keputusan")
return(hasil)
}
Ho = matrix(c(0,0,0,0),4,1)
multivariate(populasi1,populasi2,Ho)## $`Vektor rata 1`
## [1] 15.96875 15.90625 27.18750 22.75000
##
## $`matriks var-covar1`
## [,1] [,2] [,3] [,4]
## [1,] 5.192540 4.545363 6.522177 5.250000
## [2,] 4.545363 13.184476 6.760081 6.266129
## [3,] 6.522177 6.760081 28.673387 14.467742
## [4,] 5.250000 6.266129 14.467742 16.645161
##
## $`Vektor rata 2`
## [1] 12.34375 13.90625 16.71875 21.93750
##
## $`matriks var-covar2`
## [,1] [,2] [,3] [,4]
## [1,] 9.136089 7.549395 4.454637 4.151210
## [2,] 7.549395 18.603831 9.650202 5.445565
## [3,] 4.454637 9.650202 28.789315 12.917339
## [4,] 4.151210 5.445565 12.917339 27.995968
##
## $`s-pooled`
## [,1] [,2] [,3] [,4]
## [1,] 7.164315 6.047379 5.488407 4.700605
## [2,] 6.047379 15.894153 8.205141 5.855847
## [3,] 5.488407 8.205141 28.731351 13.692540
## [4,] 4.700605 5.855847 13.692540 22.320565
##
## $`T-square`
## [,1]
## [1,] 98.22791
##
## $`F tabel`
## [1] 2.527907
##
## $`c-square`
## [1] 10.62578
##
## $keputusan
## [,1]
## [1,] "Tolak Ho"#Selang Kepercayaan Simultan
msimultan<-function(mdt1, mdt2,alfa=0.05){
#populasi 1
n1<-nrow(mdt1)
p1<-ncol(mdt1)
#Membentuk vektor rata-rata sampel
xbar1=apply(mdt1, 2, mean)
varcovar1=cov(mdt1)
#populasi 2
n2<-nrow(mdt2)
p2<-ncol(mdt2)
#Membentuk vektor rata-rata sampel
xbar2=apply(mdt2, 2, mean)
varcovar2=cov(mdt2)
selisih_xbar=matrix(c(xbar1-xbar2),n1,1)
#membentuk S_pooled
spl=((n1-1)/(n1+n2-2))*varcovar1 + ((n2-1)/(n1+n2-2))*varcovar2
Ftab<-qf(1-alfa,p1,n1+n2-p1-1)
csq<-((n1+n2-2)*p1/(n1+n2-p1-1))*Ftab
#Definisi variabel loop
ci_lower=matrix()
ci_upper=matrix()
ci=matrix()
for (i in 1 : p1){
ci_lower=selisih_xbar[i,]-sqrt(csq)*sqrt(((1/n1)+(1/n2))*spl[i,i])
ci_upper=selisih_xbar[i,]+sqrt(csq)*sqrt(((1/n1)+(1/n2))*spl[i,i])
ci<-c(ci_lower, ci_upper)
print(ci)
}
}
msimultan(populasi1,populasi2)## [1] 1.443739 5.806261
## [1] -1.24892 5.24892
## [1] 6.100592 14.836908
## [1] -3.037608 4.662608#Otomatis
library(MVTests)##
## Attaching package: 'MVTests'## The following object is masked from 'package:datasets':
##
## iris#Pengujian
TwoSamplesHT2(data=gabung,group = ind_populasi, alpha = 0.05, Homogenity = TRUE)## $HT2
## [,1]
## [1,] 98.22791
##
## $F
## [,1]
## [1,] 23.36874
##
## $df
## [1] 4 59
##
## $p.value
## [,1]
## [1,] 1.307589e-11
##
## $CI
## Lower Upper Important Variables?
## x1 1.443739 5.806261 *TRUE*
## x2 -1.248920 5.248920 FALSE
## x3 6.100592 14.836908 *TRUE*
## x4 -3.037608 4.662608 FALSE
##
## $alpha
## [1] 0.05
##
## $Descriptive1
## x1 x2 x3 x4
## Means 15.968750 15.906250 27.187500 22.750000
## Sd 2.278715 3.631043 5.354754 4.079848
##
## $Descriptive2
## x1 x2 x3 x4
## Means 12.343750 13.906250 16.718750 21.937500
## Sd 3.022596 4.313216 5.365567 5.291122
##
## $Test
## [1] "TwoSamplesHT2"
##
## attr(,"class")
## [1] "MVTests" "list"#Selang Bonferroni
mbonferroni<-function(mdt1, mdt2,alfa=0.05){
#populasi 1
n1<-nrow(mdt1)
p1<-ncol(mdt1)
#Membentuk vektor rata-rata sampel
xbar1=apply(mdt1, 2, mean)
varcovar1=cov(mdt1)
#populasi 2
n2<-nrow(mdt2)
p2<-ncol(mdt2)
#Membentuk vektor rata-rata sampel
xbar2=apply(mdt2, 2, mean)
varcovar2=cov(mdt2)
selisih_xbar=matrix(c(xbar1-xbar2),n1,1)
#membentuk S_pooled
spl=((n1-1)/(n1+n2-2))*varcovar1 + ((n2-1)/(n1+n2-2))*varcovar2
t_alfa=alfa/(2*p1)
#Definisi variabel loop
ci_lower=matrix()
ci_upper=matrix()
ci=matrix()
for (i in 1 : p1){
ci_lower=selisih_xbar[i,]-(qt(1-t_alfa,n1+n2-2))*sqrt(((1/n1)+(1/n2))*spl[i,i])
ci_upper=selisih_xbar[i,]+(qt(1-t_alfa,n1+n2-2))*sqrt(((1/n1)+(1/n2))*spl[i,i])
ci<-c(ci_lower, ci_upper)
print(ci)
}
}
mbonferroni(populasi1,populasi2)## [1] 1.903487 5.346513
## [1] -0.56414 4.56414
## [1] 7.021275 13.916225
## [1] -2.226115 3.851115#Asumsi Sigma1 tidak sama dengan Sigma2
multvar<-function(mdt1, mdt2, Ho,alfa=0.05){
#populasi 1
n1<-nrow(mdt1)
p1<-ncol(mdt1)
#Membentuk vektor rata-rata sampel
xbar1=apply(mdt1, 2, mean)
varcovar1=cov(mdt1)
#populasi 2
n2<-nrow(mdt2)
p2<-ncol(mdt2)
#Membentuk vektor rata-rata sampel
xbar2=apply(mdt2, 2, mean)
varcovar2=cov(mdt2)
#T hotelling
Tsq<-t((xbar1-xbar2)-Ho)%*%solve((varcovar1/n1)+(varcovar2/n2))%*%((xbar1-xbar2)-Ho)
#wilayah kritis
chis<-qchisq(1-alfa,p1)
uji<-ifelse(Tsq>chis, "Tolak Ho","Gagal Tolak Ho")
hasil<-list(xbar1, varcovar1, xbar2, varcovar2,Tsq,chis,uji)
names(hasil)<-c("Vektor rata 1","matriks var-covar1","Vektor rata 2",
"matriks var-covar2","Statistik Uji","Chi Tabel","Keputusan")
return(hasil)
}
Ho = matrix(c(0,0,0,0),4,1)
multvar(populasi1,populasi2,Ho)## $`Vektor rata 1`
## [1] 15.96875 15.90625 27.18750 22.75000
##
## $`matriks var-covar1`
## [,1] [,2] [,3] [,4]
## [1,] 5.192540 4.545363 6.522177 5.250000
## [2,] 4.545363 13.184476 6.760081 6.266129
## [3,] 6.522177 6.760081 28.673387 14.467742
## [4,] 5.250000 6.266129 14.467742 16.645161
##
## $`Vektor rata 2`
## [1] 12.34375 13.90625 16.71875 21.93750
##
## $`matriks var-covar2`
## [,1] [,2] [,3] [,4]
## [1,] 9.136089 7.549395 4.454637 4.151210
## [2,] 7.549395 18.603831 9.650202 5.445565
## [3,] 4.454637 9.650202 28.789315 12.917339
## [4,] 4.151210 5.445565 12.917339 27.995968
##
## $`Statistik Uji`
## [,1]
## [1,] 98.22791
##
## $`Chi Tabel`
## [1] 9.487729
##
## $Keputusan
## [,1]
## [1,] "Tolak Ho"library("DescTools")
HotellingsT2Test(gabung~ind_populasi,test="chi")##
## Hotelling's two sample T2-test
##
## data: gabung by ind_populasi
## T.2 = 98.228, df = 4, p-value < 2.2e-16
## alternative hypothesis: true location difference is not equal to c(0,0,0,0)#Deteksi dan Uji Kesamaan Varians 2 Populasi
#Membangun fungsi Box'M
n1 = nrow(populasi1)
n2 = nrow(populasi2)
p = 4 #jumlah variabel
g = 2 #jumlah populasi
u = (sum(1/(n1-1),sum(1/(n2-1)))-1/sum(n1-1,n2-1))*((2*p^2+3*p-1)/(6*(p+1)*(g-1)))
spl=((n1-1)/(n1+n2-2))*var(populasi1) + ((n2-1)/(n1+n2-2))*var(populasi2)
M = sum(n1-1,n2-1)*log(det(spl))-sum((n1-1)*log(det(var(populasi1))),(n2-1)*log(det(var(populasi2))))
C = (1-u)*M
C## [1] 13.72729#Membangun fungsi boxM F
df1 = n1-1
df2 = n2-1
C2 =((p-1)*(p+2))/(6*(g-1))*(sum(1/df1^2,1/df2^2)-1/(sum(df1,df2))^2)
a1 = 1/2*(g-1)*p*(p+1)
a2 = (a1+2)/abs(C2-u^2)
b1=(1-u-a1/a2)/a1
b2 = (1-u-2/a2)/a2
lnM = 1/2*sum((n1-1)*log(det(var(populasi1))),(n2-1)*log(det(var(populasi2))))-1/2*sum(n1-1,n2-1)*log(det(spl))
#Kondisi C2 dan C1
if (C2 > u^2){
F = -2*b1*lnM
print(F)
} else {
F = -(a2*b2*lnM)/(a1*(1+2*b2*lnM))
print(F)
}## [1] 1.371926#Keputusan
qf(0.05,a1,a2)## [1] 0.3939864#Box'M
library(biotools)## Loading required package: MASS## ---
## biotools version 4.2boxM(gabung,ind_populasi)##
## Box's M-test for Homogeneity of Covariance Matrices
##
## data: gabung
## Chi-Sq (approx.) = 13.727, df = 10, p-value = 0.1858#Coba di library (MVnTests)