rm(list=ls())

library(spdep)
library(lmtest)
library(spatialreg)
library(spData)

# Input Data
data <- read.csv(file.choose(),header=T,sep=",",dec=".",fileEncoding = 'UTF-8-BOM')
data
# Variabel-variabel penelitian
Y = data$Curah.Hujan
X1 = data$Tekanan.Permukaan
X2 = data$Suhu.Udara
X3 = data$Kelembapan
X4 = data$Kecepatan.Angin
X5 = data$Penyinaran.Matahari
data2 = data.frame(Y, X1, X2, X3, X4, X5)
data2

# Statistik Deskriptif
summary(data2)
sd(data2$Y)
sd(data2$X1)
sd(data2$X2)
sd(data2$X3)
sd(data2$X4)
sd(data2$X5)

var(data2$Y)
var(data2$X1)
var(data2$X2)
var(data2$X3)
var(data2$X4)
var(data2$X5)

# Nilai Korelasi
cor(data2)

# Hitogram dan Plot QQ-Normal
hist(data2$Y)
hist(data2$X1)
hist(data2$X2)
hist(data2$X3)
hist(data2$X4)
hist(data2$X5)

qqnorm(data2$Y)
qqline(data2$Y)

qqnorm(data2$X1)
qqline(data2$X1)

qqnorm(data2$X2)
qqline(data2$X2)

qqnorm(data2$X3)
qqline(data2$X3)

qqnorm(data2$X4)
qqline(data2$X4)

qqnorm(data2$X5)
qqline(data2$X5)

# 1. Menghitung Matriks Bobot Spasial Invers Jarak 
normXY = function(X,Y){
  norm = matrix(NA, nrow = length(X), ncol = length(Y))
  for (i in 1:length(X)){
    for (j in 1:length(Y)){
      norm[i,j] = sqrt((X[i]-X[j])^2+(Y[i]-Y[j])^2)
    }
  }
  return(norm)
}
weighted = function(X,Y){
  weight = matrix(NA, nrow = length(X), ncol = length(Y))
  norm = normXY(X,Y)
  for (i in 1:length(X)){
    for (j in 1:length(Y)){
      if(i==j){
        weight[i,j] = as.numeric(0)
      }else{
        weight[i,j] = 1/(norm[i,j])
      }
    }
  }
  return(weight)
}

# Matriks Bobot Spasial Invers Jarak 
bobot_invers_jarak <- weighted(data$Lattitude..derajat.,data$Longitude..derajat.)
bobot_invers_jarak

# Matriks Bobot Spasial Invers Jarak Terstandardisasi
bobot_invers_jarak_terstandarisasi <-bobot_invers_jarak[,c(1:51)]/(rowSums(bobot_invers_jarak[,c(1:51)]))
bobot_invers_jarak_terstandarisasi

# Mengonversi matriks bobot invers jarak terstandarisasi menjadi objek listw
matbot <- mat2listw(bobot_invers_jarak_terstandarisasi, style = "W")

# 2. Ideks Moran
# Melakukan uji Moran
moran.test(data2$Y, matbot, zero.policy = TRUE)
moran.test(data2$X1, matbot, zero.policy = TRUE)
moran.test(data2$X2, matbot, zero.policy = TRUE)
moran.test(data2$X3, matbot, zero.policy = TRUE)
moran.test(data2$X4, matbot, zero.policy = TRUE)
moran.test(data2$X5, matbot, zero.policy = TRUE)

# Moran Scatterplot
moran.plot(data2$Y,matbot ,zero.policy=TRUE)
moran.plot(data2$X1,matbot,zero.policy=TRUE)
moran.plot(data2$X2,matbot,zero.policy=TRUE)
moran.plot(data2$X3,matbot,zero.policy=TRUE)
moran.plot(data2$X4,matbot,zero.policy=TRUE)
moran.plot(data2$X5,matbot,zero.policy=TRUE)

# 3. Menaksir Parameter Model Linear
# Menaksir Parameter Linear
reg1 <- lm(Y ~ X1, data = data2)
reg2 <- lm(Y ~ X2, data = data2)
reg3 <- lm(Y ~ X3, data = data2)
reg4 <- lm(Y ~ X4, data = data2)
reg5 <- lm(Y ~ X5, data = data2)

reg6 <- lm(Y ~ X1 + X2, data = data2)
reg7 <- lm(Y ~ X1 + X3, data = data2)
reg8 <- lm(Y ~ X1 + X4, data = data2)
reg9 <- lm(Y ~ X1 + X5, data = data2)
reg10 <- lm(Y ~ X2 + X3, data = data2)

reg11 <- lm(Y ~ X2 + X4, data = data2)
reg12 <- lm(Y ~ X2 + X5, data = data2)
reg13 <- lm(Y ~ X3 + X4, data = data2)
reg14 <- lm(Y ~ X3 + X5, data = data2)
reg15 <- lm(Y ~ X4 + X5, data = data2)

reg16 <- lm(Y ~ X1 + X2 + X3, data = data2)
reg17 <- lm(Y ~ X1 + X2 + X4, data = data2)
reg18 <- lm(Y ~ X1 + X2 + X5, data = data2)
reg19 <- lm(Y ~ X1 + X3 + X4, data = data2)
reg20 <- lm(Y ~ X1 + X3 + X5, data = data2)

reg21 <- lm(Y ~ X1 + X4 + X5, data = data2)
reg22 <- lm(Y ~ X2 + X3 + X4, data = data2)
reg23 <- lm(Y ~ X2 + X3 + X5, data = data2)
reg24 <- lm(Y ~ X2 + X4 + X5, data = data2)
reg25 <- lm(Y ~ X3 + X4 + X5, data = data2)

reg26 <- lm(Y ~ X1 + X2 + X3 + X4, data = data2)
reg27 <- lm(Y ~ X1 + X2 + X3 + X5, data = data2)
reg28 <- lm(Y ~ X1 + X2 + X4 + X5, data = data2)
reg29 <- lm(Y ~ X1 + X3 + X4 + X5, data = data2)
reg30 <- lm(Y ~ X2 + X3 + X4 + X5, data = data2)

reg31 <- lm(Y ~ X1 + X2 + X3 + X4 + X5, data = data2)

# Menampilkan summary dari masing-masing model regresi
summary(reg1) # pval = 0.002499
summary(reg2)  # pval = 1.01e-05
summary(reg3) # pval = 9.818e-12
summary(reg4) # pval = 0.02244
summary(reg5) # pval = 2.558e-08
summary(reg6) # pval = 9.425e-09
summary(reg7) # pval = 9.001e-12
summary(reg8) # pval = 0.00954
summary(reg9) # pval = 2.299e-10
summary(reg10) # pval = 3.147e-12
summary(reg11) # pval = 3.908e-05
summary(reg12) # pval = 5.409e-12
summary(reg13) # pval = 8.945e-12
summary(reg14) # pval = 9.788e-15
summary(reg15) # pval = 9.711e-08
summary(reg16) # pval = 2.199e-12
summary(reg17) # pval = 1.338e-08
summary(reg18) # pval = 1.452e-12
summary(reg19) # pval = 2.182e-11
summary(reg20) # pval = 7.283e-14
summary(reg21) # pval = 5.103e-10
summary(reg22) # pval = 1.16e-11
summary(reg23) # pval = 5.851e-14
summary(reg24) # pval = 9.772e-13
summary(reg25) # pval = 2.923e-15 (terkecil)
summary(reg26) # pval = 1.221e-11
summary(reg27) # pval = 1.693e-13
summary(reg28) # pval = 1.025e-13
summary(reg29) # pval = 2.13e-14
summary(reg30) # pval = 2.193e-14
summary(reg31) # pval = 1.366e-13

# 4. Menaksir Parameter Model SAR-X
lm.LMtests(reg1,matbot,test=c("LMerr","LMlag"),zero.policy=TRUE)
sar1=lagsarlm(reg1,data=data2,matbot,zero.policy=TRUE)

lm.LMtests(reg2,matbot,test=c("LMerr","LMlag"),zero.policy=TRUE)
sar2=lagsarlm(reg2,data=data2,matbot,zero.policy=TRUE)

lm.LMtests(reg3,matbot,test=c("LMerr","LMlag"),zero.policy=TRUE)
sar3=lagsarlm(reg3,data=data2,matbot,zero.policy=TRUE)

lm.LMtests(reg4,matbot,test=c("LMerr","LMlag"),zero.policy=TRUE)
sar4=lagsarlm(reg4,data=data2,matbot,zero.policy=TRUE)

lm.LMtests(reg5,matbot,test=c("LMerr","LMlag"),zero.policy=TRUE)
sar5=lagsarlm(reg5,data=data2,matbot,zero.policy=TRUE)

lm.LMtests(reg6,matbot,test=c("LMerr","LMlag"),zero.policy=TRUE)
sar6=lagsarlm(reg6,data=data2,matbot,zero.policy=TRUE)

lm.LMtests(reg7,matbot,test=c("LMerr","LMlag"),zero.policy=TRUE)
sar7=lagsarlm(reg7,data=data2,matbot,zero.policy=TRUE)

lm.LMtests(reg8,matbot,test=c("LMerr","LMlag"),zero.policy=TRUE)
sar8=lagsarlm(reg8,data=data2,matbot,zero.policy=TRUE)

lm.LMtests(reg9,matbot,test=c("LMerr","LMlag"),zero.policy=TRUE)
sar9=lagsarlm(reg9,data=data2,matbot,zero.policy=TRUE)

lm.LMtests(reg10,matbot,test=c("LMerr","LMlag"),zero.policy=TRUE)
sar10=lagsarlm(reg10,data=data2,matbot,zero.policy=TRUE)

lm.LMtests(reg11,matbot,test=c("LMerr","LMlag"),zero.policy=TRUE)
sar11=lagsarlm(reg11,data=data2,matbot,zero.policy=TRUE)

lm.LMtests(reg12,matbot,test=c("LMerr","LMlag"),zero.policy=TRUE)
sar12=lagsarlm(reg12,data=data2,matbot,zero.policy=TRUE)

lm.LMtests(reg13,matbot,test=c("LMerr","LMlag"),zero.policy=TRUE)
sar13=lagsarlm(reg13,data=data2,matbot,zero.policy=TRUE)

lm.LMtests(reg14,matbot,test=c("LMerr","LMlag"),zero.policy=TRUE)
sar14=lagsarlm(reg14,data=data2,matbot,zero.policy=TRUE)

lm.LMtests(reg15,matbot,test=c("LMerr","LMlag"),zero.policy=TRUE)
sar14=lagsarlm(reg15,data=data2,matbot,zero.policy=TRUE)

lm.LMtests(reg15, matbot, test = c("LMerr", "LMlag"), zero.policy = TRUE)
sar15 <- lagsarlm(reg15, data = data2, matbot, zero.policy = TRUE)

lm.LMtests(reg16, matbot, test = c("LMerr", "LMlag"), zero.policy = TRUE)
sar16 <- lagsarlm(reg16, data = data2, matbot, zero.policy = TRUE)

lm.LMtests(reg17, matbot, test = c("LMerr", "LMlag"), zero.policy = TRUE)
sar17 <- lagsarlm(reg17, data = data2, matbot, zero.policy = TRUE)

lm.LMtests(reg18, matbot, test = c("LMerr", "LMlag"), zero.policy = TRUE)
sar18 <- lagsarlm(reg18, data = data2, matbot, zero.policy = TRUE)

lm.LMtests(reg19, matbot, test = c("LMerr", "LMlag"), zero.policy = TRUE)
sar19 <- lagsarlm(reg19, data = data2, matbot, zero.policy = TRUE)

lm.LMtests(reg20, matbot, test = c("LMerr", "LMlag"), zero.policy = TRUE)
sar20 <- lagsarlm(reg20, data = data2, matbot, zero.policy = TRUE)

lm.LMtests(reg21, matbot, test = c("LMerr", "LMlag"), zero.policy = TRUE)
sar21 <- lagsarlm(reg21, data = data2, matbot, zero.policy = TRUE)

lm.LMtests(reg22, matbot, test = c("LMerr", "LMlag"), zero.policy = TRUE)
sar22 <- lagsarlm(reg22, data = data2, matbot, zero.policy = TRUE)

lm.LMtests(reg23, matbot, test = c("LMerr", "LMlag"), zero.policy = TRUE)
sar23 <- lagsarlm(reg23, data = data2, matbot, zero.policy = TRUE)

lm.LMtests(reg24, matbot, test = c("LMerr", "LMlag"), zero.policy = TRUE)
sar24 <- lagsarlm(reg24, data = data2, matbot, zero.policy = TRUE)

lm.LMtests(reg25, matbot, test = c("LMerr", "LMlag"), zero.policy = TRUE)
sar25 <- lagsarlm(reg25, data = data2, matbot, zero.policy = TRUE)

lm.LMtests(reg26, matbot, test = c("LMerr", "LMlag"), zero.policy = TRUE)
sar26 <- lagsarlm(reg26, data = data2, matbot, zero.policy = TRUE)

lm.LMtests(reg27, matbot, test = c("LMerr", "LMlag"), zero.policy = TRUE)
sar27 <- lagsarlm(reg27, data = data2, matbot, zero.policy = TRUE)

lm.LMtests(reg28, matbot, test = c("LMerr", "LMlag"), zero.policy = TRUE)
sar28 <- lagsarlm(reg28, data = data2, matbot, zero.policy = TRUE)

lm.LMtests(reg29, matbot, test = c("LMerr", "LMlag"), zero.policy = TRUE)
sar29 <- lagsarlm(reg29, data = data2, matbot, zero.policy = TRUE)

lm.LMtests(reg30, matbot, test = c("LMerr", "LMlag"), zero.policy = TRUE)
sar30 <- lagsarlm(reg30, data = data2, matbot, zero.policy = TRUE)

lm.LMtests(reg31, matbot, test = c("LMerr", "LMlag"), zero.policy = TRUE)
sar31 <- lagsarlm(reg31, data = data2, matbot, zero.policy = TRUE)

# Menampilkan summary dari masing-masing model SAR-X
summary(sar1)  #9.5057e-09
summary(sar2)  #1.148e-08
summary(sar3)  #8.9107e-07
summary(sar4)  #4.1507e-08
summary(sar5)  #0.00036853
summary(sar6)  #1.6696e-07
summary(sar7)  #6.366e-06
summary(sar8)  #1.0951e-08
summary(sar9)  #0.00071452
summary(sar10) #1.4625e-05
summary(sar11) #6.0351e-09 (Terkecil)
summary(sar12) #0.001055
summary(sar13) #1.1709e-06
summary(sar14) #0.00034012
summary(sar15) #0.00034012
summary(sar16) #0.0001078
summary(sar17) #7.1774e-08
summary(sar18) #0.0015122
summary(sar19) #3.8921e-06
summary(sar20) #0.00455
summary(sar21) #0.0010324
summary(sar22) #7.4225e-06
summary(sar23) #0.0051608
summary(sar24) #0.0027535
summary(sar25) #0.012105
summary(sar26) #4.7347e-05
summary(sar27) #0.010029
summary(sar28) #0.0049487
summary(sar29) #0.012297
summary(sar30) #0.011861
summary(sar31) #0.013259

# 5. Prediksi dengan SAR-X model 11 dan Mencari MAPE
# Nilai Y Topi
Y_hat = sar11$fitted.values

# Menghitung nilai absolut dari selisih antara nilai observasi dan prediksi
abs_error <- abs(data2$Y - Y_hat)

# Menghitung MAPE
MAPE <- mean(abs_error / data2$Y) * 100
MAPE

# 6. Menentukan data training (data tersampel) dan testing (data tak tersampel) untuk kriging
data_kriging <- data.frame(data$Kota.Kab, data$Lattitude, data$Longitude, Y_hat)
# Menghitung nilai rata-rata latitude dan longitude
mean_latitude <- mean(data_kriging$data.Lattitude)
mean_longitude <- mean(data_kriging$data.Longitude)

# Membuat fungsi untuk menghitung jarak Euclidean
euclidean_distance <- function(lat1, lon1, lat2, lon2) {
  sqrt((lat1 - lat2)^2 + (lon1 - lon2)^2)
}

# Hitung jarak dari setiap titik data ke pusat data
data_kriging$DistanceToCenter <- euclidean_distance(data_kriging$data.Lattitude, data_kriging$data.Longitude, mean_latitude, mean_longitude)

# Urutkan data berdasarkan jarak ke pusat data
sorted_data <- data_kriging[order(data_kriging$DistanceToCenter), ]

# Pilih 10 data pertama sebagai data testing
data_testing <- sorted_data[1:10, ]
data_testing

# Pilih data training dengan menghilangkan data testing
data_training <- sorted_data[-(1:10), ]
data_training
# Data testing :
# Kota Semarang, Kota Magelang, Kab. Kendal, Kab. Sleman, Kab. Wonosobo, Kota Yogyakarta, Kab. Kudus, Kab. Sragen, Kab. Kulonprogo, Kab. Karanganyar

# 7. Melakukan Ordinary Point Kriging
library(scatterplot3d)
library(sp)
library(gstat)

# Uji Normalitas dan Boxplot
boxplot(data_training$Y_hat)
shapiro.test(data_training$Y_hat)

# koordinat xy dan presipitasi
scatterplot3d(data_training$data.Lattitude, data_training$data.Longitude, data_training$Y_hat, pch = 16, color = "blue",
              scale.y = 1, angle = 70)

# 7.1 Semivariogram
data_kriging = data.frame(curah_hujan = data_training$Y_hat, Lat = data_training$data.Lattitude, Long = data_training$data.Longitude)
coordinates(data_kriging) = ~Lat+Long
variogram_eks = variogram(curah_hujan~1, data = data_kriging, cloud = F)
plot(variogram_eks)

# semivariogram teoretis
sill = var(data_kriging$curah_hujan)
# model spherical
sph.model = vgm(psill = sill, model = "Sph", width = 0.01)
sph.fit = fit.variogram(object = variogram_eks, model = sph.model)
sph.fit
plot(variogram_eks, pl = F, model = sph.fit, col = "blue", cex = 1, lwd = 0.5, lty = 1, pch = 20,
     main = "Variogram models - Spherical", xlab = "Distance (m)", ylab = "Semivariance")
# model eksponential
exp.model = vgm(psill = sill, model = "Exp", width = 0.001)
exp.fit = fit.variogram(object = variogram_eks, model = exp.model)
exp.fit
plot(variogram_eks, pl = F, model = exp.fit, col = "blue", cex = 1, lwd = 0.5, lty = 1, pch = 20,
     main = "Variogram models - Exponential", xlab = "Distance (m)", ylab = "Semivariance")
# model gaussian
gau.model = vgm(psill = sill, model = "Gau", width = 0.1)
gau.fit = fit.variogram(object = variogram_eks, model = gau.model)
gau.fit
plot(variogram_eks, pl = F, model = sph.fit, col = "blue", cex = 1, lwd = 0.5, lty = 1, pch = 20,
     main = "Variogram models - Gaussian", xlab = "Distance (m)", ylab = "Semivariance")

attributes(sph.fit)$SSErr
attributes(exp.fit)$SSErr
attributes(gau.fit)$SSErr

# point kriging
data_testing = data.frame(curah_hujan1 = data_testing$Y_hat, Lat1 = data_testing$data.Lattitude, Long1 = data_testing$data.Longitude)
coordinates(data_testing) = ~Lat1+Long1
interpolasi = krige(curah_hujan~1, data_kriging, data_testing, model = exp.fit)

data_aktual = data_testing$curah_hujan1
data_prediksi = interpolasi$var1.pred
hasil_prediksi = data.frame(data_aktual, data_prediksi)
hasil_prediksi

mape = mean(abs((data_aktual - data_prediksi)/data_aktual))*100
mape

rmse = sqrt(mean((data_aktual - data_prediksi)^2))
rmse