##Membaca data CSV yang akan diolah
setwd("C:/Users/USER/Documents/Kuliah/Semester 6/Analisis Model Prediktif")
suppressPackageStartupMessages(library(tidyverse))
library(readr)
AMP_Contoh <- read_delim("data gabungan.csv",
                         delim = ",", escape_double = FALSE, trim_ws = TRUE)
View(AMP_Contoh)

##standarisasi data
# Simpan nama variabel
nama_variabel <- names(AMP_Contoh)
nama_variabel

index_provinsi <- which(names(AMP_Contoh) == "Provinsi")
index_tahun <- which(names(AMP_Contoh) == "Tahun")
index_pengeluaran <- which(names(AMP_Contoh) == "Pengeluaran konsumsi")
# Buat subset data tanpa variabel "Provinsi" dan "Tahun"
AMP_Contoh_subset <- AMP_Contoh[, -c(index_provinsi, index_tahun, index_pengeluaran)] 

AMP_Contoh_numeric <- as.data.frame(lapply(AMP_Contoh_subset, as.numeric))
AMP_Contoh_standardized <- as.data.frame(scale(AMP_Contoh_numeric))

AMP_Contoh_standardized$Provinsi <- AMP_Contoh$Provinsi
AMP_Contoh_standardized$Tahun <- AMP_Contoh$Tahun
index_pengeluaran2 <- which(names(AMP_Contoh_standardized) == "Pengeluaran konsumsi")
attach(AMP_Contoh_standardized)
library(PerformanceAnalytics)
chart.Correlation(AMP_Contoh_standardized[,1:9]) #yang dipakai dari tabel 1 sampe 10

#uji bartlett
uji_bart <- function(x)
{
  method <- "Bartlett's test of sphericity"
  data.name <- deparse(substitute(x))
  x <- subset(x, complete.cases(x)) 
  n <- nrow(x)
  p <- ncol(x)
  chisq <- (1-n+(2*p+5)/6)*log(det(cor(x)))
  df <- p*(p-1)/2
  p.value <- pchisq(chisq, df, lower.tail=FALSE)
  names(chisq) <- "Khi-squared"
  names(df) <- "df"
  return(structure(list(statistic=chisq, parameter=df, p.value=p.value,
                        method=method, data.name=data.name), class="htest"))
}
uji_bart(AMP_Contoh_standardized[,1:9])



##Uji KMO
library(psych)
kmo <- function(x)
{
  x <- subset(x, complete.cases(x))       # menghilangkan data kosong (NA)
  r <- cor(x)                             # Membuat matrix korelasi
  r2 <- r^2                               # nilai koefisien untuk r squared
  i <- solve(r)                           # Inverse matrix dari matrix korelasi
  d <- diag(i)                            # element diagonal dari inverse matrix
  p2 <- (-i/sqrt(outer(d, d)))^2          # koefisien korelasi Parsial kuadrat
  diag(r2) <- diag(p2) <- 0               # menghapus element diagonal 
  KMO <- sum(r2)/(sum(r2)+sum(p2))
  MSA <- colSums(r2)/(colSums(r2)+colSums(p2))
  return(list(KMO=KMO, MSA=MSA))
}
kmo(AMP_Contoh_standardized[,1:9]) #sampel cukup untuk digunakan jika kmo>0.5

###Analisis PCA###
library(factoextra)
pca_AMP <- prcomp(AMP_Contoh_standardized[,1:9], scale. = TRUE)

pca_AMP
round(pca_AMP$rotation)

#bayaknya jumlah faktor dilihat dari nilai eigen value > 1 atau nilai varians > 80% 
round(pca_AMP$sdev^2,2) 
fviz_eig(pca_AMP, addlabels = TRUE, ylim = c(0, 80))
PCA_AMP_FIX <- round(pca_AMP$rotation[,1:4],2)
PCA_AMP_FIX
variabel_PCA<- as.data.frame(PCA_AMP_FIX)
variabel_PCA
#akses sanitasi, rata2 luas hunian, kepemilikan bpjs non iuran, jumlah kendaraan bermotor


#Membuat Variabel Biner
median_kesejahteraan <- median(AMP_Contoh$`Pengeluaran konsumsi`)
median_kesejahteraan
AMP_Contoh$sejahtera_binary <- as.factor(ifelse(AMP_Contoh$`Pengeluaran konsumsi` > median_kesejahteraan, 1, 0))


# Memisahkan variabel yang diminta
akses_sanitasi <- AMP_Contoh_standardized$Akses.Sanitasi
sumber_penerangan <- AMP_Contoh_standardized$Sumber.Penerangan
kepemilikan_bpjs <- AMP_Contoh_standardized$Kepemilikan.bpjs.non
jumlah_kendaraan <- AMP_Contoh_standardized$Jumlah.kendaran.bermotor
Kesejahteraan <- AMP_Contoh$sejahtera_binary

# Membuat data frame baru dengan variabel-variabel yang dipisahkan
data_baru <- data.frame(Kesejahteraan = Kesejahteraan,
                        Akses.Sanitasi = akses_sanitasi,
                        Sumber_penerangan = sumber_penerangan,
                        Kepemilikan_BPJS = kepemilikan_bpjs,
                        Jumlah_Kendaraan = jumlah_kendaraan)

#Memeriksa data frame baru
head(data_baru)

#--------------------knn--------------------
library(caTools)
library(dplyr)
library(ggplot2)
library(caret)
library(class)
library(readr)

#memilih tipe data pada masing" variabel
sapply(data_baru, class)
#dimensi dataset
dim(data_baru)
#Melakukan standarisasi data 
standard.features <- scale(data_baru[,2:5])
#Split data > data train & data test
set.seed(200)
dat.samp<-sample.split(data_baru$Kesejahteraan, SplitRatio = 0.70)
dat.train<-subset(data_baru, dat.samp==TRUE)
dat.test<-subset(data_baru, dat.samp==FALSE)
dim(dat.train)
dat.samp

#membuat prediksi dengan menggunakan knn 
#k=2
predicted.type <-  knn(dat.train[2:5],dat.test[2:5],dat.train$Kesejahteraan,k=2)
#Error di prediksi
error <- mean(predicted.type!=dat.test$Kesejahteraan)
error
#Confusion Matrix
confusionMatrix(predicted.type,dat.test$Kesejahteraan)

#Mengecek nilai error dari masing-masing k
predicted.type <- NULL
error.rate <- NULL
for (i in 1:15) {
  predicted.type <- knn(dat.train[2:5],dat.test[2:5],dat.train$Kesejahteraan,k=i)
  error.rate[i] <- mean(predicted.type!=dat.test$Kesejahteraan)
  
}
knn.error <- as.data.frame(cbind(k=1:15,error.type =error.rate))
ggplot(knn.error,aes(k,error.type))+ 
  geom_point()+ 
  geom_line() + 
  scale_x_continuous(breaks=1:15)+ 
  theme_bw() +
  xlab("Value of K") +
  ylab('Error')

#k=14
predicted.type <-  knn(dat.train[2:5],dat.test[2:5],dat.train$Kesejahteraan,k=14)
#Error di prediksi
error <- mean(predicted.type!=dat.test$Kesejahteraan)
error
#Confusion Matrix
confusionMatrix(predicted.type,dat.test$Kesejahteraan)

#k=6
predicted.type <-  knn(dat.train[2:5],dat.test[2:5],dat.train$Kesejahteraan,k=6)
#Error di prediksi
error <- mean(predicted.type!=dat.test$Kesejahteraan)
error
#Confusion Matrix
confusionMatrix(predicted.type,dat.test$Kesejahteraan)


#--------------------random forest--------------------
library(randomForest)
library(datasets)
library(caret)

str(data_baru)
#membagi data
set.seed(222)
#Split data > data train & data test
set.seed(200)
dat.samp<-sample.split(data_baru$Kesejahteraan, SplitRatio = 0.70)
dat.train<-subset(data_baru, dat.samp==TRUE)
dat.test<-subset(data_baru, dat.samp==FALSE)
dim(dat.train)
dim(dat.test)
dat.samp

rf <- randomForest(Kesejahteraan~., data=dat.train, proximity=TRUE) 
print(rf)

p1 <- predict(rf, dat.train)
confusionMatrix(p1, dat.train$Kesejahteraan)

p2 <- predict(rf, dat.test)
confusionMatrix(p2, dat.test$Kesejahteraan)


#--------------------naive bayes--------------------
#installing required packages
# importing the required libraries
library(mlbench)
library(caret)
library(e1071)


#Split data > data train & data test
set.seed(200)
dat.samp<-sample.split(data_baru$Kesejahteraan, SplitRatio = 0.70)
dat.train<-subset(data_baru, dat.samp==TRUE)
dat.test<-subset(data_baru, dat.samp==FALSE)
dim(dat.train)
dat.samp

# training the model
model <- naiveBayes(Kesejahteraan ~ ., data = dat.train)
# making predictions
predictions <- predict(model, newdata = dat.test)

# evaluating the model
confusion_matrix <- table(predictions, dat.test$Kesejahteraan)
accuracy <- sum(diag(confusion_matrix)) / sum(confusion_matrix)
print(confusion_matrix)
cat("Accuracy:", round(accuracy * 100, 2), "%\n")