Langkah 1 : Import Data

library(readxl)
## Warning: package 'readxl' was built under R version 4.4.3
heart_dataset <- read_excel("C:/Users/denih/OneDrive/Dokumen/1. Statmul/3. TUGAS FIX/Tugas 7/heart_dataset.xlsx")
head(heart_dataset)

Langkah 2: Memilih Variabel untuk Analisis

data_numeric <- heart_dataset[, sapply(heart_dataset, is.numeric)]  
data_numeric
X1 <- heart_dataset$age
X2 <- heart_dataset$trestbps
X3 <- heart_dataset$chol
Y1 <- heart_dataset$thalach
Y2 <- heart_dataset$oldpeak

data1 <- data.frame(X1,X2,X3,Y1,Y2)
data1
X <- data1[, c("X1", "X2", "X3")] 
Y <- data1[, c("Y1", "Y2")]

# Standarisasi X
X_scaled <- scale(X)

# Standarisasi Y
Y_scaled <- scale(Y)

# Menggabungkan kembali dalam satu data frame
data_standarisasi <- data.frame(X_scaled, Y_scaled)

Langkah 3: Melihat 10 Data Pertama

head(data_standarisasi, 10)

Langkah 4: Pemeriksaan Data

  1. Melihat Tipe Data
sapply(data_standarisasi, class)
##        X1        X2        X3        Y1        Y2 
## "numeric" "numeric" "numeric" "numeric" "numeric"
  1. Cek Missing Value
colSums(is.na(data_standarisasi))
## X1 X2 X3 Y1 Y2 
##  0  0  0  0  0
  1. Melihat Dimensi Data
dim(data_standarisasi)
## [1] 303   5

Langkah 5:

a) Menghitung Rata-Rata

matriks_data <- as.matrix(data_standarisasi,303,5)
x_bar <- colMeans(matriks_data)
x_bar
##            X1            X2            X3            Y1            Y2 
##  7.841179e-17 -6.870192e-16 -9.761112e-17 -5.879854e-16 -3.476317e-17

b) Membuat Matriks Varian-Kovarian

cov_matriks <- cov(matriks_data)
cov_matriks
##            X1          X2           X3           Y1          Y2
## X1  1.0000000  0.27935091  0.213677957 -0.398521938  0.21001257
## X2  0.2793509  1.00000000  0.123174207 -0.046697728  0.19321647
## X3  0.2136780  0.12317421  1.000000000 -0.009939839  0.05395192
## Y1 -0.3985219 -0.04669773 -0.009939839  1.000000000 -0.34418695
## Y2  0.2100126  0.19321647  0.053951920 -0.344186948  1.00000000

Langkah 6: Uji Asumsi

a) Uji Linearitas

pairs(data_standarisasi,
      panel = function(x, y) {
        points(x, y, pch = 16, col = "blue")  
        abline(lm(y ~ x), col = "red") 
      },
      main = "Scatter Plot Matrix")

# b) Uji Normalitas Multivariat

Di <- mahalanobis(matriks_data,x_bar,cov_matriks)
Di
##   [1]  2.6619901 13.5605334  3.3472857  2.9999638  5.6356193  2.1043348
##   [7]  1.0472588  2.3714771  7.9136130  6.2987876  0.5474786  2.6453894
##  [13]  1.1383230  5.1724071  1.8645751  1.4025256  5.7989989  4.4715757
##  [19]  4.2143730  3.8546563  1.1576717  2.1614226  3.5067953  1.6905583
##  [25]  4.8433524  7.2703356  2.4084574  4.9514102 11.5502124  0.9762128
##  [31]  4.0538635  5.3174559  3.3588133  0.6411601  2.1628518  3.3617050
##  [37]  2.5843470  2.1149120  2.8933388  7.2563899  2.5931720  1.8707538
##  [43]  7.8099222  0.8174728  7.6727667  4.0730132  3.0963586  2.0190370
##  [49]  5.2370274  1.1853853  0.4996580  4.2458807  1.5286331  6.5202015
##  [55]  3.6470114  0.9657868  3.0017894  3.4464728  5.8381237  2.2865633
##  [61]  7.9128869  4.2615110  5.6518298  6.5751942  2.5275463  7.3978660
##  [67]  3.6142941  1.6705252  2.0032087  3.7130410  0.8319066  4.9743887
##  [73]  9.9645701  3.1315610  1.9551719  0.5732583  3.2026894  2.3677970
##  [79]  3.3556472  3.3231935  3.5935341  3.6979289  7.1021165  3.9813110
##  [85]  8.0694869 41.7742020  4.5561005  3.9577870  2.9841485  5.1826634
##  [91]  1.7143931  2.4133969  1.6440204  1.5520778  5.2891480  7.1137011
##  [97]  8.4252534  2.2958100  5.2950388  1.4862227  4.1811787 12.7691369
## [103]  6.4811233  4.7237426  1.7606959  5.1397768  6.8453765  2.2780806
## [109]  0.8618069  2.0926260 11.3642938  9.6077515  3.9321335  2.8189368
## [115]  0.9640409  3.9518440  4.0332570  3.6154161  3.6147578  2.3378195
## [121]  3.1308597  3.7476838  3.6804359  3.1871592  6.7885151  6.5760649
## [127]  2.7210747  6.0081922  2.2018173  7.9333270  5.2904424  1.5141895
## [133]  4.0083370  3.5262076  4.9172754  1.3451324 10.4386169  2.9351361
## [139]  3.4983073  6.3061625  1.7127172  5.4756511  2.6363304  6.5694018
## [145]  8.4633372  6.1413039  2.1032899  2.8851862  1.9208564  4.0439265
## [151]  4.3548792 11.2506557  6.5288534  3.6463849  4.9945630  2.1349337
## [157]  2.0337847  5.1623178  1.1275289  1.4809772  2.0492187  4.0991824
## [163]  5.4180595  5.6883032  5.6883032  6.3382979  5.0305252  7.4066571
## [169]  1.2233041  4.9001415  0.4379082  2.3580797  2.6933628  7.8766782
## [175]  2.5644913 10.0178631  2.8049805  4.9920751  7.4928905  5.9484354
## [181]  5.4004580  4.5293557  4.7736393  6.2090327  4.4356148  4.0041093
## [187]  0.5216601  4.5242520  0.9090685  4.3817587  2.0783182  1.8241484
## [193]  4.3469532  2.9048730  6.1240953  9.3761491  7.9487618  4.6200470
## [199]  5.7114683  5.1151920  3.5883067  3.1657794  5.4200516  8.4304042
## [205] 26.4909387  0.8934729  2.3974314  3.3621052  2.8327710  4.1996162
## [211]  0.6691993  6.8423344  4.3886282  1.9589694  0.2748986 10.7731750
## [217]  5.8987579  3.6342747  3.2003094  1.5930577 16.6149623 16.8082308
## [223]  4.6098907 20.6370322  5.0113367  7.3356597  5.3880976  8.3181993
## [229]  6.4981574  3.1310105  2.5272067  6.4368242  3.8050784  6.3433974
## [235]  6.5406289  2.9710605  3.1514626  2.5340408 11.5098562  7.1017630
## [241]  6.3007288  8.3191753  2.5916082 11.2315317  5.5775371  1.2902209
## [247] 10.6406020  7.4653749 17.4403523  3.4245151 10.1424445  3.5158976
## [253]  4.5017677  8.2335259  6.7218882  5.7050449  3.2765606  3.8833793
## [259]  1.8284992 14.5177901  9.4544070  1.8715225  8.1958452  7.3861807
## [265]  7.4504322  5.6883791 14.8576926  5.4201430  6.2499916  5.4196525
## [271]  1.5761917  2.5268871 13.9300540  4.7990695  5.9593145  1.3423289
## [277]  5.0425405  1.0956580  3.0725718  8.0623669  8.7894713  0.7651637
## [283]  1.9051378  6.2481680  1.8502407  7.4813094  1.9548472  3.3589614
## [289]  8.8117599  1.8363807  3.9836767 13.3802319  6.6961055  4.2882373
## [295]  6.4946791  9.6906951  3.9585102 14.5973770  3.3069792  4.0862172
## [301]  8.1306945  7.4374511  2.2824397
hasil <- data.frame(Obs = 1:length(Di), Mahalanobis_Distance = Di)

hasil<- hasil[order(hasil$Mahalanobis_Distance), ]
hasil$Rank <- c(1:303)

hasil$Probability<-((hasil$Rank-0.5)/303)
hasil
hasil$X2<-qchisq(hasil$Probability,5)
hasil
plot(hasil$Mahalanobis_Distance,hasil$X2,
     xlab= "Mahalanobis Distance",ylab="Chi-Square",,col="red",main = "QQ plot Normalitas",
     pch=19,
     cex=0.8)

c) Uji Non-Multikolineritas

VIF <- function(x){
  VIF <- diag(solve(cor(x)))
  result <- ifelse(VIF > 10,"multicolinearity","non
                   multicolinearity")
  data1 <- data.frame(VIF,result)
  return(data1)
}
VIF(data_standarisasi)

Langkah 7: Membuat Partisi 2 Himpunan Matriks Kovarian

P11 <- cov_matriks[1:3, 1:3]  # Kovarians antara variabel X (3 variabel)
P11
##           X1        X2        X3
## X1 1.0000000 0.2793509 0.2136780
## X2 0.2793509 1.0000000 0.1231742
## X3 0.2136780 0.1231742 1.0000000
P12 <- cov_matriks[1:3, 4:5]  # Kovarians antara variabel X dan Y
P12
##              Y1         Y2
## X1 -0.398521938 0.21001257
## X2 -0.046697728 0.19321647
## X3 -0.009939839 0.05395192
P21 <- cov_matriks[4:5, 1:3]  # Kovarians antara variabel Y dan X
P21
##            X1          X2           X3
## Y1 -0.3985219 -0.04669773 -0.009939839
## Y2  0.2100126  0.19321647  0.053951920
P22 <- cov_matriks[4:5, 4:5]  # Kovarians antara variabel Y (2 variabel)
P22
##            Y1         Y2
## Y1  1.0000000 -0.3441869
## Y2 -0.3441869  1.0000000

Langkah 8: Mencari Invers Matriks Akar Kuadrat

a) Mencari Nilai Σ11^-1/2

eig.P11 <- eigen(P11)
nilai.eigen.P11 <- eig.P11$values
nilai.eigen.P11
## [1] 1.4172628 0.8823956 0.7003416
l1.11 <- nilai.eigen.P11[1]
l2.11 <- nilai.eigen.P11[2]
l3.11 <- nilai.eigen.P11[3]
vektor.eigen.P11 <- eig.P11$vectors
vektor.eigen.P11
##            [,1]       [,2]       [,3]
## [1,] -0.6438029 -0.1173564  0.7561384
## [2,] -0.5787731 -0.5717129 -0.5815205
## [3,] -0.5005393  0.8120171 -0.3001477
v1.11<-matrix(vektor.eigen.P11[,1]) 
v1.11
##            [,1]
## [1,] -0.6438029
## [2,] -0.5787731
## [3,] -0.5005393
v2.11<-matrix(vektor.eigen.P11[,2]) 
v2.11
##            [,1]
## [1,] -0.1173564
## [2,] -0.5717129
## [3,]  0.8120171
v3.11<-matrix(vektor.eigen.P11[,3]) 
v3.11
##            [,1]
## [1,]  0.7561384
## [2,] -0.5815205
## [3,] -0.3001477
sig11 <- ((v1.11 %*% t(v1.11)) / sqrt(l1.11)) +
         ((v2.11 %*% t(v2.11)) / sqrt(l2.11)) +
         ((v3.11 %*% t(v3.11)) / sqrt(l3.11))
sig11
##            [,1]        [,2]        [,3]
## [1,]  1.0460227 -0.14100622 -0.10195602
## [2,] -0.1410062  1.03342203 -0.04229889
## [3,] -0.1019560 -0.04229889  1.02003923

b) Mencari Nilai Σ22^-1/2

eig.P22 <- eigen(P22)
nilai.eigen.P22 <- eig.P22$values
nilai.eigen.P22
## [1] 1.3441869 0.6558131
l1.22 <- nilai.eigen.P22[1]
l2.22 <- nilai.eigen.P22[2]
vektor.eigen.P22 <- eig.P22$vectors
vektor.eigen.P22
##            [,1]       [,2]
## [1,] -0.7071068 -0.7071068
## [2,]  0.7071068 -0.7071068
v1.22 <- matrix(vektor.eigen.P22[,1])
v1.22
##            [,1]
## [1,] -0.7071068
## [2,]  0.7071068
v2.22 <- matrix(vektor.eigen.P22[,2])
v2.22
##            [,1]
## [1,] -0.7071068
## [2,] -0.7071068
sig22 <- ((v1.22 %*% t(v1.22)) / sqrt(l1.22)) +
         ((v2.22 %*% t(v2.22)) / sqrt(l2.22))
sig22
##          [,1]     [,2]
## [1,] 1.048680 0.186158
## [2,] 0.186158 1.048680

Langkah 9: Menghitung Matriks M dan N

a) Untuk Σ11

M.sig.11 <- sig11 %*% P12 %*% solve(P22) %*% P21 %*% sig11
M.sig.11
##             [,1]        [,2]         [,3]
## [1,]  0.16990531 0.005498550 -0.011374069
## [2,]  0.00549855 0.033106656  0.007361137
## [3,] -0.01137407 0.007361137  0.002575673
eigen11 <- eigen(M.sig.11)
nilai.eigenM11 <- eigen11$values
nilai.eigenM11
## [1] 1.708561e-01 3.473157e-02 2.387979e-18
l1.M11<-nilai.eigenM11[1]
l1.M11
## [1] 0.1708561
l2.M11<-nilai.eigenM11[2]
l2.M11
## [1] 0.03473157
l3.M11<-nilai.eigenM11[3]
l3.M11
## [1] 2.387979e-18
vektor.eigen.M11<-eigen11$vectors
vektor.eigen.M11
##             [,1]        [,2]        [,3]
## [1,]  0.99717204 -0.02023499  0.07237727
## [2,]  0.03628729  0.97300437 -0.22791605
## [3,] -0.06581152  0.22989789  0.97098703

Mencari nilai & vektor eigen a untuk Σ11

e1<-matrix(vektor.eigen.M11[,1])
e1
##             [,1]
## [1,]  0.99717204
## [2,]  0.03628729
## [3,] -0.06581152
e2<-matrix(vektor.eigen.M11[,2])
e2
##             [,1]
## [1,] -0.02023499
## [2,]  0.97300437
## [3,]  0.22989789
e3<-matrix(vektor.eigen.M11[,3])
e3
##             [,1]
## [1,]  0.07237727
## [2,] -0.22791605
## [3,]  0.97098703

b) Untuk Σ22

M.sig.22 <- sig22%*%P21%*%solve(P11)%*%P12%*%sig22
M.sig.22
##             [,1]        [,2]
## [1,]  0.15864667 -0.03889639
## [2,] -0.03889639  0.04694097
eigen22 <- eigen(M.sig.22)
nilai.eigenM22 <- eigen22$values
nilai.eigenM22
## [1] 0.17085607 0.03473157
l1.M22<-nilai.eigenM22[1]
l1.M22
## [1] 0.1708561
l2.M22<-nilai.eigenM22[2]
l2.M22
## [1] 0.03473157
vektor.eigen.M22<-eigen22$vectors
vektor.eigen.M22
##            [,1]       [,2]
## [1,] -0.9541002 -0.2994877
## [2,]  0.2994877 -0.9541002

Mencari nilai & vektor eigen a untuk Σ22

v1.M22<-matrix(vektor.eigen.M22[,1])
v1.M22
##            [,1]
## [1,] -0.9541002
## [2,]  0.2994877
v2.M22<-matrix(vektor.eigen.M22[,2])
v2.M22
##            [,1]
## [1,] -0.2994877
## [2,] -0.9541002

Langkah 10: Mencari Koefisien A dan B

a) Koefisien A

a1<-sig11%*%e1
a1
##            [,1]
## [1,]  1.0446577
## [2,] -0.1003236
## [3,] -0.1703329
a2<-sig11%*%e2
a2
##            [,1]
## [1,] -0.1818054
## [2,]  0.9986530
## [3,]  0.1954109
a3<-sig11%*%e3
a3
##              [,1]
## [1,]  0.008847877
## [2,] -0.286810786
## [3,]  0.992706153

b) Koefisien B

b1<-sig22%*%v1.M22
b1
##            [,1]
## [1,] -0.9447937
## [2,]  0.1364534
b2<-sig22%*%v2.M22
b2
##            [,1]
## [1,] -0.4916801
## [2,] -1.0562978

Langkah 11: Menghitung Korelasi Pasangan Kanonik

# Korelasi pasangan kanonik pertama
U1V1 <- (t(a1) %*% P12 %*% b1) / ((sqrt(t(a1) %*% P11 %*% a1)) * (sqrt(t(b1) %*% P22 %*% b1)))
U1V1
##           [,1]
## [1,] 0.4133474
# Korelasi pasangan kanonik kedua
U2V2 <- (t(a2) %*% P12 %*% b2) / ((sqrt(t(a2) %*% P11 %*% a2)) * (sqrt(t(b2) %*% P22 %*% b2)))
U2V2
##            [,1]
## [1,] -0.1863641
CC1 <- cancor(X[, 1:3], Y[, 1:2])  # Hanya menggunakan 3 variabel dari X dan 2 dari Y
CC1
## $cor
## [1] 0.4133474 0.1863641
## 
## $xcoef
##             [,1]          [,2]          [,3]
## X1  0.0066188758 -0.0011519059  5.605951e-05
## X2 -0.0003291669  0.0032766310 -9.410407e-04
## X3 -0.0001891070  0.0002169491  1.102122e-03
## 
## $ycoef
##            [,1]        [,2]
## Y1 -0.002373560 0.001235224
## Y2  0.006762706 0.052350713
## 
## $xcenter
##        X1        X2        X3 
##  54.36634 131.62376 246.26403 
## 
## $ycenter
##         Y1         Y2 
## 149.646865   1.039604

Memunculkan korelasi kanonik

canonical.correlation <- CC1$cor
canonical.correlation 
## [1] 0.4133474 0.1863641

memunculkan raw cannonical

raw.cannonical <- CC1$ycoef
raw.cannonical
##            [,1]        [,2]
## Y1 -0.002373560 0.001235224
## Y2  0.006762706 0.052350713
raw.cannonical <- CC1$xcoef
raw.cannonical
##             [,1]          [,2]          [,3]
## X1  0.0066188758 -0.0011519059  5.605951e-05
## X2 -0.0003291669  0.0032766310 -9.410407e-04
## X3 -0.0001891070  0.0002169491  1.102122e-03

Langkah 12: Uji Signifikansi Korelasi Kanonik

a) Uji Secara Simultan

library(stats)

simultan.cor <- cancor(X_scaled, Y_scaled)

# Menampilkan hasil korelasi kanonik secara simultan
cat("Korelasi kanonik secara simultan:\n")
## Korelasi kanonik secara simultan:
print(simultan.cor$cor)
## [1] 0.4133474 0.1863641
# Uji signifikansi secara simultan (Wilks' Lambda Test)
n <- nrow(X_scaled)  # Jumlah sampel
p <- ncol(X_scaled)  # Jumlah variabel X
q <- ncol(Y_scaled)  # Jumlah variabel Y
m <- min(p, q)       # Jumlah pasangan korelasi kanonik

# Menghitung Wilks' Lambda
wilks_lambda <- prod(1 - simultan.cor$cor^2)
df1 <- (p * q)
df2 <- (n - 1) - (p + q + 1) / 2
F_stat <- ((1 - wilks_lambda) / wilks_lambda) * (df2 / df1)
p_value <- pf(F_stat, df1, df2, lower.tail = FALSE)

cat("\nWilks' Lambda:", wilks_lambda, "\n")
## 
## Wilks' Lambda: 0.8003465
cat("F-statistic:", F_stat, "\n")
## F-statistic: 12.43137
cat("p-value:", p_value, "\n")
## p-value: 1.66211e-12
if (p_value < 0.05) {
  cat("Kesimpulan: Tolak H0, ada paling tidak satu korelasi kanonik yang signifikan.\n")
} else {
  cat("Kesimpulan: Terima H0, tidak ada korelasi kanonik yang signifikan.\n")
}
## Kesimpulan: Tolak H0, ada paling tidak satu korelasi kanonik yang signifikan.

b) Uji Secara Parsial

parsial.cor_Y1 <- cancor(X_scaled, Y_scaled[, 1, drop = FALSE])
print("Korelasi kanonik secara parsial antara Y1 dan X1-X3:")
## [1] "Korelasi kanonik secara parsial antara Y1 dan X1-X3:"
print(parsial.cor_Y1$cor)
## [1] 0.4106358
parsial.cor_Y2 <- cancor(X_scaled, Y_scaled[, 2, drop = FALSE])
print("Korelasi kanonik secara parsial antara Y2 dan X1-X3:")
## [1] "Korelasi kanonik secara parsial antara Y2 dan X1-X3:"
print(parsial.cor_Y2$cor)
## [1] 0.2524701

Langkah 13: Interpretasi Fungsi Kanonik

a) Menghitung Bobot Kanonik

bobot_X <- CC1$xcoef[, 1, drop = FALSE]  # Ambil pasangan kanonik pertama
bobot_Y <- CC1$ycoef[, 1, drop = FALSE]

print("Bobot Kanonik untuk Variabel X :")
## [1] "Bobot Kanonik untuk Variabel X :"
print(bobot_X)
##             [,1]
## X1  0.0066188758
## X2 -0.0003291669
## X3 -0.0001891070
print("Bobot Kanonik untuk Variabel Y:")
## [1] "Bobot Kanonik untuk Variabel Y:"
print(bobot_Y)
##            [,1]
## Y1 -0.002373560
## Y2  0.006762706

b) Menghitung Muatan Kanonik

bobot_X <- as.matrix(bobot_X)
bobot_Y <- as.matrix(bobot_Y)

muatan_X <- cor(X_scaled, X_scaled %*% bobot_X)
muatan_Y <- cor(Y_scaled, Y_scaled %*% bobot_Y)


print("Muatan Kanonik untuk Variabel X:")
## [1] "Muatan Kanonik untuk Variabel X:"
print(muatan_X)
##         [,1]
## X1 0.9983176
## X2 0.2303257
## X3 0.1823264
print("Muatan Kanonik untuk Variabel Y:")
## [1] "Muatan Kanonik untuk Variabel Y:"
print(muatan_Y)
##         [,1]
## Y1 -0.595052
## Y2  0.959392

c) Menentukan Variabel Paling Berpengaruh

max_X <- which.max(abs(muatan_X[, 1]))  # Pasangan kanonik pertama
max_Y <- which.max(abs(muatan_Y[, 1]))

print("Variabel X yang paling berpengaruh dalam pasangan kanonik pertama:")
## [1] "Variabel X yang paling berpengaruh dalam pasangan kanonik pertama:"
print(colnames(X)[max_X])
## [1] "X1"
print("Variabel Y yang paling berpengaruh dalam pasangan kanonik pertama:")
## [1] "Variabel Y yang paling berpengaruh dalam pasangan kanonik pertama:"
print(colnames(Y)[max_Y])
## [1] "Y2"