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)
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)
head(data_standarisasi, 10)
sapply(data_standarisasi, class)
## X1 X2 X3 Y1 Y2
## "numeric" "numeric" "numeric" "numeric" "numeric"
colSums(is.na(data_standarisasi))
## X1 X2 X3 Y1 Y2
## 0 0 0 0 0
dim(data_standarisasi)
## [1] 303 5
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
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
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)
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)
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
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
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
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
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
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
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
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
b1<-sig22%*%v1.M22
b1
## [,1]
## [1,] -0.9447937
## [2,] 0.1364534
b2<-sig22%*%v2.M22
b2
## [,1]
## [1,] -0.4916801
## [2,] -1.0562978
# 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
canonical.correlation <- CC1$cor
canonical.correlation
## [1] 0.4133474 0.1863641
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
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.
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
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
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
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"