Import Library
library(readxl)
library(psych)
## Warning: package 'psych' was built under R version 4.5.2
library(GPArotation)
##
## Attaching package: 'GPArotation'
## The following objects are masked from 'package:psych':
##
## equamax, varimin
library(corrplot)
## corrplot 0.95 loaded
library(knitr)
library(kableExtra)
## Warning: package 'kableExtra' was built under R version 4.5.2
library(dplyr)
##
## Attaching package: 'dplyr'
## The following object is masked from 'package:kableExtra':
##
## group_rows
## The following objects are masked from 'package:stats':
##
## filter, lag
## The following objects are masked from 'package:base':
##
## intersect, setdiff, setequal, union
Load Dataset
data1 <- read_excel("G:/Other computers/My Laptop/Tugas/Semester 5/StatMul/project/data1.xlsx")
data2 <- read_excel("G:/Other computers/My Laptop/Tugas/Semester 5/StatMul/project/data2.xlsx")
head(data1)
## # A tibble: 6 × 12
## No `Kabupaten/Kota` X1 X2 X3 X4 X5 X6 X7 X8 X9
## <dbl> <chr> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
## 1 1 Nias 0 55 36.0 38 0 32.9 7.7 90.7 13
## 2 2 Mandailing Natal 2 131 23.0 0 1 58.3 30.7 70.0 140
## 3 3 Tapanuli Selatan 1 0 18.8 12 0 52.4 24.4 62.4 152
## 4 4 Tapanuli Tengah 4 92 23.0 22 0 62.4 22.7 71.0 41
## 5 5 Tapanuli Utara 1 89 14.2 62 0 61.2 19.2 57.9 22
## 6 6 Toba Samosir 3 88 21.8 64 1 67.6 14.8 53.8 20
## # ℹ 1 more variable: X10 <dbl>
head(data2)
## # A tibble: 6 × 12
## No `Kabupaten/Kota` X1 X2 X3 X4 X5 X6 X7 X8 X9
## <dbl> <chr> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
## 1 1 Nias 58.8 32.2 33.5 5.12 89.5 89.4 74.3 49.5 63.2
## 2 2 Mandailing Natal 46.7 9.4 8.48 1.35 55.1 72.8 53.7 10.3 62.5
## 3 3 Tapanuli Selatan 42.8 10.1 5.61 1.63 61.1 72.7 61.6 7.69 62.9
## 4 4 Tapanuli Tengah 54.9 6.22 6.31 1.68 46.0 71.8 59.1 9.52 59.7
## 5 5 Tapanuli Utara 37.3 2.58 1.2 0.79 26.3 75.4 72.9 7.79 62.9
## 6 6 Toba Samosir 27.5 1.81 3.33 1.66 16.7 65.6 59.2 5.02 55.1
## # ℹ 1 more variable: X10 <dbl>
Select Numeric Variables Only
data1_num <- data1 %>% select(where(is.numeric))
data2_num <- data2 %>% select(where(is.numeric))
Merge Dataset
data <- rbind(data1_num, data2_num)
dim(data)
## [1] 66 11
Descriptive Statistics
psych::describe(data)
## vars n mean sd median trimmed mad min max range
## No 1 66 17.00 9.59 17.00 17.00 11.86 1.00 33.00 32.00
## X1 2 66 21.69 19.01 22.28 20.24 27.10 0.00 64.00 64.00
## X2 3 66 75.77 123.98 22.76 49.46 32.67 0.00 721.00 721.00
## X3 4 66 14.88 12.12 19.27 14.51 16.77 0.00 35.97 35.97
## X4 5 66 66.95 159.62 5.25 28.71 5.46 0.00 997.00 997.00
## X5 6 66 15.78 24.49 5.00 10.52 6.60 0.00 90.09 90.09
## X6 7 66 59.70 15.69 59.36 59.95 18.09 26.16 89.44 63.28
## X7 8 66 32.91 19.31 25.55 31.38 15.34 6.20 75.90 69.70
## X8 9 66 37.62 27.48 39.72 36.77 38.50 0.44 90.70 90.26
## X9 10 66 58.69 43.04 56.89 52.16 11.22 4.00 265.00 261.00
## X10 11 66 89628.05 142844.17 157.50 60395.11 230.54 0.00 738440.00 738440.00
## skew kurtosis se
## No 0.00 -1.26 1.18
## X1 0.39 -1.20 2.34
## X2 2.76 9.64 15.26
## X3 -0.01 -1.62 1.49
## X4 3.90 17.25 19.65
## X5 1.78 2.08 3.01
## X6 -0.14 -0.91 1.93
## X7 0.72 -0.68 2.38
## X8 0.11 -1.37 3.38
## X9 2.22 7.10 5.30
## X10 2.44 6.94 17582.90
Correlation Matrix Visualization
corrplot(cor(data), method = "color", tl.cex = 0.8)
KMO Test of Sampling Adequacy
KMO(data)
## Kaiser-Meyer-Olkin factor adequacy
## Call: KMO(r = data)
## Overall MSA = 0.71
## MSA for each item =
## No X1 X2 X3 X4 X5 X6 X7 X8 X9 X10
## 0.40 0.74 0.62 0.69 0.59 0.57 0.86 0.84 0.75 0.36 0.80
Bartlett’s Test of Sphericity
cortest.bartlett(cor(data), n = nrow(data))
## $chisq
## [1] 524.0892
##
## $p.value
## [1] 1.021485e-77
##
## $df
## [1] 55
Parallel Analysis to Determine Number of Factors
fa.parallel(data, fa="fa")
## Warning in fa.stats(r = r, f = f, phi = phi, n.obs = n.obs, np.obs = np.obs, :
## The estimated weights for the factor scores are probably incorrect. Try a
## different factor score estimation method.
## Warning in fa.stats(r = r, f = f, phi = phi, n.obs = n.obs, np.obs = np.obs, :
## The estimated weights for the factor scores are probably incorrect. Try a
## different factor score estimation method.
## Warning in fa.stats(r = r, f = f, phi = phi, n.obs = n.obs, np.obs = np.obs, :
## The estimated weights for the factor scores are probably incorrect. Try a
## different factor score estimation method.
## Warning in fa.stats(r = r, f = f, phi = phi, n.obs = n.obs, np.obs = np.obs, :
## The estimated weights for the factor scores are probably incorrect. Try a
## different factor score estimation method.
## Warning in fac(r = r, nfactors = nfactors, n.obs = n.obs, rotate = rotate, : An
## ultra-Heywood case was detected. Examine the results carefully
## Parallel analysis suggests that the number of factors = 3 and the number of components = NA
Factor Analysis
fa_result <- fa(data, nfactors = 4, rotate = "varimax", fm = "pa")
## maximum iteration exceeded
## Warning in fa.stats(r = r, f = f, phi = phi, n.obs = n.obs, np.obs = np.obs, :
## The estimated weights for the factor scores are probably incorrect. Try a
## different factor score estimation method.
## Warning in fac(r = r, nfactors = nfactors, n.obs = n.obs, rotate = rotate, : An
## ultra-Heywood case was detected. Examine the results carefully
fa_result
## Factor Analysis using method = pa
## Call: fa(r = data, nfactors = 4, rotate = "varimax", fm = "pa")
## Standardized loadings (pattern matrix) based upon correlation matrix
## PA1 PA3 PA2 PA4 h2 u2 com
## No -0.05 -0.02 0.05 -0.56 0.32 0.682 1.0
## X1 0.81 -0.52 -0.02 -0.06 0.93 0.072 1.7
## X2 -0.10 0.30 0.75 -0.10 0.68 0.325 1.4
## X3 -0.08 0.94 0.24 0.06 0.95 0.047 1.2
## X4 -0.10 0.12 1.10 0.14 1.25 -0.250 1.1
## X5 0.98 0.06 0.03 0.08 0.97 0.027 1.0
## X6 0.61 -0.45 -0.21 0.34 0.73 0.270 2.8
## X7 0.73 -0.24 -0.13 0.20 0.64 0.361 1.4
## X8 -0.22 0.89 0.04 0.14 0.87 0.134 1.2
## X9 0.09 0.00 0.21 0.28 0.13 0.873 2.1
## X10 0.19 -0.54 -0.19 0.21 0.41 0.594 1.9
##
## PA1 PA3 PA2 PA4
## SS loadings 2.63 2.61 1.97 0.65
## Proportion Var 0.24 0.24 0.18 0.06
## Cumulative Var 0.24 0.48 0.66 0.72
## Proportion Explained 0.33 0.33 0.25 0.08
## Cumulative Proportion 0.33 0.67 0.92 1.00
##
## Mean item complexity = 1.5
## Test of the hypothesis that 4 factors are sufficient.
##
## df null model = 55 with the objective function = 8.66 with Chi Square = 524.09
## df of the model are 17 and the objective function was 0.49
##
## The root mean square of the residuals (RMSR) is 0.03
## The df corrected root mean square of the residuals is 0.05
##
## The harmonic n.obs is 66 with the empirical chi square 5.72 with prob < 0.99
## The total n.obs was 66 with Likelihood Chi Square = 28.39 with prob < 0.041
##
## Tucker Lewis Index of factoring reliability = 0.917
## RMSEA index = 0.1 and the 90 % confidence intervals are 0.021 0.165
## BIC = -42.84
## Fit based upon off diagonal values = 0.99
Factor Loadings Table
loadings <- as.data.frame(fa_result$loadings[1:ncol(data), ])
kable(round(loadings, 3), caption = "Factor Loadings After Varimax Rotation") %>%
kable_styling(full_width = FALSE)
Factor Loadings After Varimax Rotation
|
|
PA1
|
PA3
|
PA2
|
PA4
|
|
No
|
-0.054
|
-0.021
|
0.050
|
-0.559
|
|
X1
|
0.810
|
-0.518
|
-0.018
|
-0.060
|
|
X2
|
-0.102
|
0.301
|
0.751
|
-0.101
|
|
X3
|
-0.084
|
0.941
|
0.236
|
0.057
|
|
X4
|
-0.098
|
0.115
|
1.098
|
0.144
|
|
X5
|
0.981
|
0.059
|
0.029
|
0.076
|
|
X6
|
0.608
|
-0.452
|
-0.208
|
0.336
|
|
X7
|
0.726
|
-0.241
|
-0.129
|
0.196
|
|
X8
|
-0.216
|
0.893
|
0.039
|
0.142
|
|
X9
|
0.092
|
-0.003
|
0.205
|
0.277
|
|
X10
|
0.187
|
-0.538
|
-0.191
|
0.212
|
Factor Diagram
fa.diagram(fa_result)
Factor Scores
factor_scores <- as.data.frame(fa_result$scores)
head(factor_scores)
## PA1 PA3 PA2 PA4
## 1 -0.7232539 1.61541328 -0.21556989 0.2737752
## 2 -0.6388670 1.09155352 -1.99555604 0.5855036
## 3 -0.7074717 0.59430467 -0.85599095 0.9644460
## 4 -0.6693256 0.73010369 -0.85580107 0.7348257
## 5 -0.7336325 0.02818056 -0.11907745 0.7977766
## 6 -0.6403298 0.35206975 -0.02135831 0.8164102