Untitled

Analisis faktor,PCA,dan Bioplot

library(readxl)
library(psych)
library(GPArotation)

## 
## Attaching package: 'GPArotation'

## The following objects are masked from 'package:psych':
## 
##     equamax, varimin

library(nFactors)
library(factoextra)

## Loading required package: ggplot2

## 
## Attaching package: 'ggplot2'

## The following objects are masked from 'package:psych':
## 
##     %+%, alpha

## Welcome! Want to learn more? See two factoextra-related books at https://goo.gl/ve3WBa

library(FactoMineR)
library(dplyr)

## 
## Attaching package: 'dplyr'

## The following objects are masked from 'package:stats':
## 
##     filter, lag

## The following objects are masked from 'package:base':
## 
##     intersect, setdiff, setequal, union

Memanggil data yang di perlukan

data <- read_excel("C:/Users/rafif faheem/Downloads/Nilai Impor Migas Bulanan Menurut Negara Asal, 2025.xlsx")

## New names:
## • `` -> `...2`
## • `` -> `...3`
## • `` -> `...4`
## • `` -> `...5`
## • `` -> `...6`
## • `` -> `...7`
## • `` -> `...8`
## • `` -> `...9`
## • `` -> `...10`
## • `` -> `...11`
## • `` -> `...12`
## • `` -> `...13`
## • `` -> `...14`

head(data)

## # A tibble: 6 × 14
##   Negara Asal Impor Miga…¹ ...2  ...3  ...4  ...5  ...6  ...7  ...8  ...9  ...10
##   <chr>                    <chr> <chr> <chr> <chr> <chr> <chr> <chr> <chr> <chr>
## 1 <NA>                     Nila… <NA>  <NA>  <NA>  <NA>  <NA>  <NA>  <NA>  <NA> 
## 2 <NA>                     2025  <NA>  <NA>  <NA>  <NA>  <NA>  <NA>  <NA>  <NA> 
## 3 <NA>                     Janu… Febr… Maret April Mei   Juni  Juli  Agus… Sept…
## 4 ASEAN :                  1384… 1378… 1336… 1277… 1086… 1134… 9917… -     -    
## 5 a. Thailand              5579  7745  7330… 4617… 7260… 8126… 9778… -     -    
## 6 b. Singapura             8136… 8758… 9554… 9443… 7596… 7325… 6062… -     -    
## # ℹ abbreviated name: ¹`Negara Asal Impor Migas`
## # ℹ 4 more variables: ...11 <chr>, ...12 <chr>, ...13 <chr>, ...14 <chr>

str(data)

## tibble [35 × 14] (S3: tbl_df/tbl/data.frame)
##  $ Negara Asal Impor Migas: chr [1:35] NA NA NA "ASEAN :" ...
##  $ ...2                   : chr [1:35] "Nilai Impor Migas Bulanan Menurut Negara Asal (Ribu US$)" "2025" "Januari" "1384649" ...
##  $ ...3                   : chr [1:35] NA NA "Februari" "1378071" ...
##  $ ...4                   : chr [1:35] NA NA "Maret" "1336632.2" ...
##  $ ...5                   : chr [1:35] NA NA "April" "1277605.2" ...
##  $ ...6                   : chr [1:35] NA NA "Mei" "1086889.7" ...
##  $ ...7                   : chr [1:35] NA NA "Juni" "1134266.4" ...
##  $ ...8                   : chr [1:35] NA NA "Juli" "991713.6" ...
##  $ ...9                   : chr [1:35] NA NA "Agustus" "-" ...
##  $ ...10                  : chr [1:35] NA NA "September" "-" ...
##  $ ...11                  : chr [1:35] NA NA "Oktober" "-" ...
##  $ ...12                  : chr [1:35] NA NA "November" "-" ...
##  $ ...13                  : chr [1:35] NA NA "Desember" "-" ...
##  $ ...14                  : chr [1:35] NA NA "Tahunan" "8589827.2" ...

A.Analisis Faktor

1.Cleaning data dan jadikan data Numerik

colnames(data)[1] <- "Negara"

# 4. Buang baris agregat (misalnya "ASEAN :", dst.)
data <- subset(data, !grepl(":", Negara))

# 5. Ubah tanda "-" jadi NA
data[data == "-"] <- NA

# 6. Ubah semua kolom (selain Negara) ke numeric
for(i in 2:ncol(data)){
  data[[i]] <- as.numeric(data[[i]])
}

## Warning: NAs introduced by coercion
## Warning: NAs introduced by coercion
## Warning: NAs introduced by coercion
## Warning: NAs introduced by coercion
## Warning: NAs introduced by coercion
## Warning: NAs introduced by coercion
## Warning: NAs introduced by coercion
## Warning: NAs introduced by coercion
## Warning: NAs introduced by coercion
## Warning: NAs introduced by coercion
## Warning: NAs introduced by coercion
## Warning: NAs introduced by coercion
## Warning: NAs introduced by coercion

# 7. Hapus kolom bulan yang seluruh nilainya nol
data <- data %>%
  select(-Negara) %>%                                # buang kolom Negara
  select(where(is.numeric)) %>%                      # pastikan hanya numerik
  select(where(~ sum(., na.rm = TRUE) > 0))          # buang kolom nol semua

# 8. Reset rownames (hilangkan nama negara yang mungkin nyangkut)
rownames(data) <- NULL

data[is.na(data)] <- 0


# 9. Cek hasil
str(data)

## tibble [33 × 8] (S3: tbl_df/tbl/data.frame)
##  $ ...2 : num [1:33] 0 2025 0 5579 813654 ...
##  $ ...3 : num [1:33] 0 0 0 7745 875850 ...
##  $ ...4 : num [1:33] 0 0 0 7330 955425 ...
##  $ ...5 : num [1:33] 0 0 0 4617 944349 ...
##  $ ...6 : num [1:33] 0 0 0 7260 759606 ...
##  $ ...7 : num [1:33] 0 0 0 8127 732559 ...
##  $ ...8 : num [1:33] 0 0 0 9779 606270 ...
##  $ ...14: num [1:33] 0 0 0 50437 5687712 ...

head(data)

## # A tibble: 6 × 8
##     ...2   ...3     ...4    ...5    ...6    ...7    ...8     ...14
##    <dbl>  <dbl>    <dbl>   <dbl>   <dbl>   <dbl>   <dbl>     <dbl>
## 1      0      0      0        0       0       0       0        0  
## 2   2025      0      0        0       0       0       0        0  
## 3      0      0      0        0       0       0       0        0  
## 4   5579   7745   7330.    4617.   7260.   8127.   9779.   50437. 
## 5 813654 875850 955425.  944349. 759606. 732559. 606270. 5687712  
## 6      0      0      0.1      0       0       0       0        0.1

2.Uji KMO

KMO(data)

## Kaiser-Meyer-Olkin factor adequacy
## Call: KMO(r = data)
## Overall MSA =  0.5
## MSA for each item = 
##  ...2  ...3  ...4  ...5  ...6  ...7  ...8 ...14 
##  0.50  0.50  0.50  0.49  0.50  0.50  0.49  0.50

3.Uji Bartlett

cortest.bartlett(data)

## R was not square, finding R from data

## $chisq
## [1] 1472.856
## 
## $p.value
## [1] 4.565659e-293
## 
## $df
## [1] 28

4.Uji Parallel

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 fac(r = r, nfactors = nfactors, n.obs = n.obs, rotate = rotate, : An
## ultra-Heywood case was detected.  Examine the results carefully

## 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

## 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

## 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

## 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

## 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

## 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

## 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

## 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

## 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.

## Parallel analysis suggests that the number of factors =  1  and the number of components =  NA

5.Uji faktor EFA

fa_result <- fa(data, nfactors = 2, rotate = "varimax")
print(fa_result, cut = 0.3)

## Factor Analysis using method =  minres
## Call: fa(r = data, nfactors = 2, rotate = "varimax")
## Standardized loadings (pattern matrix) based upon correlation matrix
##        MR1  MR2   h2      u2 com
## ...2  0.78 0.61 0.99 0.00604 1.9
## ...3  0.75 0.65 1.00 0.00185 2.0
## ...4  0.75 0.66 1.00 0.00348 2.0
## ...5  0.75 0.65 0.98 0.01757 2.0
## ...6  0.75 0.66 0.99 0.00599 2.0
## ...7  0.75 0.66 0.99 0.00688 2.0
## ...8  0.66 0.75 1.00 0.00490 2.0
## ...14 0.75 0.66 1.00 0.00037 2.0
## 
##                        MR1  MR2
## SS loadings           4.42 3.53
## Proportion Var        0.55 0.44
## Cumulative Var        0.55 0.99
## Proportion Explained  0.56 0.44
## Cumulative Proportion 0.56 1.00
## 
## Mean item complexity =  2
## Test of the hypothesis that 2 factors are sufficient.
## 
## df null model =  28  with the objective function =  51.68 with Chi Square =  1472.86
## df of  the model are 13  and the objective function was  16.05 
## 
## The root mean square of the residuals (RMSR) is  0 
## The df corrected root mean square of the residuals is  0 
## 
## The harmonic n.obs is  33 with the empirical chi square  0.01  with prob <  1 
## The total n.obs was  33  with Likelihood Chi Square =  435.9  with prob <  5.7e-85 
## 
## Tucker Lewis Index of factoring reliability =  0.338
## RMSEA index =  0.992  and the 90 % confidence intervals are  0.928 1.091
## BIC =  390.44
## Fit based upon off diagonal values = 1
## Measures of factor score adequacy             
##                                                    MR1  MR2
## Correlation of (regression) scores with factors   0.96 0.94
## Multiple R square of scores with factors          0.92 0.89
## Minimum correlation of possible factor scores     0.84 0.79

6.Diagram Faktor Analisis

fa.diagram(fa_result)

B.PCA

1.Summary Data(Ringkasan Data)

summary(data)

##       ...2              ...3              ...4                ...5          
##  Min.   :      0   Min.   :      0   Min.   :      0.0   Min.   :      0.0  
##  1st Qu.:      5   1st Qu.:      0   1st Qu.:      0.0   1st Qu.:      9.8  
##  Median :    407   Median :     89   Median :    276.3   Median :    444.9  
##  Mean   : 150530   Mean   : 173784   Mean   : 189548.9   Mean   : 152698.0  
##  3rd Qu.:  27548   3rd Qu.:  47847   3rd Qu.:  49248.3   3rd Qu.:  57436.9  
##  Max.   :2482727   Max.   :2867440   Max.   :3127556.7   Max.   :2519517.9  
##       ...6                ...7                ...8          
##  Min.   :      0.0   Min.   :      0.0   Min.   :      0.0  
##  1st Qu.:      0.6   1st Qu.:      0.1   1st Qu.:      0.0  
##  Median :    201.5   Median :    302.4   Median :    220.9  
##  Mean   : 160222.5   Mean   : 134652.3   Mean   : 152331.9  
##  3rd Qu.:  26002.7   3rd Qu.:  26897.6   3rd Qu.:   9778.9  
##  Max.   :2643671.6   Max.   :2221763.4   Max.   :2513476.6  
##      ...14          
##  Min.   :0.000e+00  
##  1st Qu.:1.440e+01  
##  Median :2.826e+03  
##  Mean   :1.114e+06  
##  3rd Qu.:2.952e+05  
##  Max.   :1.838e+07

Menjalankan Fungsi PCA

# reset rownames
rownames(data) <- NULL

# 3. Jalankan PCA (standardisasi sangat penting)
pca_result <- prcomp(data, scale. = TRUE)

3.Ringkasan dan Loading Matrix PCA

summary(pca_result)

## Importance of components:
##                           PC1    PC2     PC3     PC4     PC5     PC6     PC7
## Standard deviation     2.8179 0.1574 0.13981 0.10888 0.05186 0.01623 0.01078
## Proportion of Variance 0.9926 0.0031 0.00244 0.00148 0.00034 0.00003 0.00001
## Cumulative Proportion  0.9926 0.9957 0.99813 0.99962 0.99995 0.99999 1.00000
##                              PC8
## Standard deviation     9.786e-05
## Proportion of Variance 0.000e+00
## Cumulative Proportion  1.000e+00

pca_result$rotation

##             PC1          PC2          PC3         PC4         PC5         PC6
## ...2  0.3531863  0.422552128  0.416085055 -0.29673273 -0.46393745  0.27596085
## ...3  0.3545691  0.060485490  0.265826497 -0.06875653 -0.15815295 -0.30631952
## ...4  0.3542899  0.077705714 -0.018444367  0.50095070 -0.03764284 -0.65767510
## ...5  0.3523319  0.233034767 -0.803189129  0.11156239 -0.21763524  0.28231005
## ...6  0.3539126 -0.028402627  0.297943606  0.49484860  0.49495041  0.52288177
## ...7  0.3537741  0.093952095 -0.146498676 -0.61120954  0.63633321 -0.18828545
## ...8  0.3514811 -0.864761576 -0.021885630 -0.15734728 -0.24692449  0.09274196
## ...14 0.3548687  0.001239247  0.006730027  0.02478007 -0.01000173 -0.01790611
##                PC7        PC8
## ...2   0.357544542  0.1276778
## ...3  -0.809264883  0.1465517
## ...4   0.397467863  0.1587906
## ...5  -0.142496740  0.1289522
## ...6  -0.061497696  0.1346492
## ...7   0.133268276  0.1133406
## ...8   0.131912487  0.1282092
## ...14 -0.004787715 -0.9343253

4.Scree Plot

fviz_eig(pca_result)

## Warning in geom_bar(stat = "identity", fill = barfill, color = barcolor, :
## Ignoring empty aesthetic: `width`.

5.Diagram Variabel PCA

fviz_pca_var(pca_result, repel = TRUE,
             col.var = "contrib") + 
  theme_minimal()

## Warning: Using `size` aesthetic for lines was deprecated in ggplot2 3.4.0.
## ℹ Please use `linewidth` instead.
## ℹ The deprecated feature was likely used in the ggpubr package.
##   Please report the issue at <https://github.com/kassambara/ggpubr/issues>.
## This warning is displayed once every 8 hours.
## Call `lifecycle::last_lifecycle_warnings()` to see where this warning was
## generated.

## Warning: `aes_string()` was deprecated in ggplot2 3.0.0.
## ℹ Please use tidy evaluation idioms with `aes()`.
## ℹ See also `vignette("ggplot2-in-packages")` for more information.
## ℹ The deprecated feature was likely used in the factoextra package.
##   Please report the issue at <https://github.com/kassambara/factoextra/issues>.
## This warning is displayed once every 8 hours.
## Call `lifecycle::last_lifecycle_warnings()` to see where this warning was
## generated.

6.Individu PCA

fviz_pca_ind(pca_result, repel = TRUE,
             geom.ind = "point",
             col.ind = "cos2") +
  theme_minimal()

C.Biplot

1.Biplot I

fviz_pca_biplot(pca_result, repel = TRUE,
                col.var = "blue", col.ind = "gray40") +
  theme_minimal()

2.Biplot II

# Jalankan PCA (dengan data yang sudah dibersihkan)
pca_result <- prcomp(data, scale. = TRUE)

# Biplot: gabungkan variabel (bulan) & observasi (negara/baris data)
fviz_pca_biplot(
  pca_result,
  repel = TRUE,              # biar label tidak tumpang tindih
  col.var = "blue",          # warna variabel (bulan)
  col.ind = "gray40",        # warna observasi (negara)
  geom.ind = "point",        # observasi jadi titik
  label = "all"              # tampilkan label semua
) + theme_minimal()

Untitled

2025-09-24