library(bayesQR)
## Warning: package 'bayesQR' was built under R version 4.5.2
dataku<-read.csv("D:/DATA ANALISIS.csv")
head(dataku)
## Kabupaten.Kota Y X1 X2 X3 X4 X5 X6
## 1 Pacitan 106 1.41 1.56 13.08 4.34 86.62 71.49
## 2 Ponorogo 271 2.30 4.19 9.11 4.74 78.75 73.70
## 3 Trenggalek 162 1.78 3.90 10.50 4.71 80.08 72.47
## 4 Tulungagung 511 2.66 4.12 6.28 4.86 75.57 75.13
## 5 Blitar 276 3.02 4.77 8.16 4.44 72.36 73.44
## 6 Kediri 616 4.04 5.10 9.95 4.95 71.31 75.18
summary(dataku)
## Kabupaten.Kota Y X1 X2
## Length:38 Min. : 106.0 Min. :0.330 Min. :1.560
## Class :character 1st Qu.: 328.0 1st Qu.:1.705 1st Qu.:3.368
## Mode :character Median : 522.5 Median :2.650 Median :4.075
## Mean : 808.0 Mean :2.631 Mean :4.041
## 3rd Qu.: 851.0 3rd Qu.:3.237 3rd Qu.:4.735
## Max. :6381.0 Max. :6.990 Max. :6.490
## X3 X4 X5 X6
## Min. : 3.060 Min. :1.670 Min. :67.52 Min. :66.72
## 1st Qu.: 6.518 1st Qu.:4.650 1st Qu.:70.97 1st Qu.:71.69
## Median : 9.215 Median :4.840 Median :73.36 Median :74.56
## Mean : 9.782 Mean :4.631 Mean :73.67 Mean :75.31
## 3rd Qu.:12.043 3rd Qu.:5.150 3rd Qu.:75.35 3rd Qu.:78.62
## Max. :20.830 Max. :5.760 Max. :86.62 Max. :84.69
#Model BayesQR (tau=0,1)
model_bqr1<-bayesQR(formula=Y~X1+X2+X3+X4+X5+X6, data=dataku, quantile=0.1, ndraw=5000)
## Current iteration :
## [1] 500
## Current iteration :
## [1] 1000
## Current iteration :
## [1] 1500
## Current iteration :
## [1] 2000
## Current iteration :
## [1] 2500
## Current iteration :
## [1] 3000
## Current iteration :
## [1] 3500
## Current iteration :
## [1] 4000
## Current iteration :
## [1] 4500
## Current iteration :
## [1] 5000
summary(model_bqr1)
##
## Type of dependent variable: continuous
## Lasso variable selection: no
## Normal approximation of posterior: yes
## Estimated quantile: 0.1
## Lower credible bound: 0.025
## Upper credible bound: 0.975
## Number of burnin draws: 0
## Number of retained draws: 5000
##
##
## Summary of the estimated beta:
##
## Bayes Estimate lower upper adj.lower adj.upper
## (Intercept) -2.89 -22.91 16.303 -41.75 36.0
## X1 115.81 100.67 128.824 -856.85 1088.5
## X2 -26.41 -33.33 -19.134 -323.10 270.3
## X3 -4.57 -8.30 -0.424 -206.32 197.2
## X4 -9.59 -16.29 -1.794 -305.00 285.8
## X5 -8.43 -9.48 -7.350 -32.33 15.5
## X6 11.42 10.50 12.523 -5.11 27.9
#Model BayesQR (tau=0,25)
model_bqr2<-bayesQR(formula=Y~X1+X2+X3+X4+X5+X6, data=dataku, quantile=0.25, ndraw=5000)
## Current iteration :
## [1] 500
## Current iteration :
## [1] 1000
## Current iteration :
## [1] 1500
## Current iteration :
## [1] 2000
## Current iteration :
## [1] 2500
## Current iteration :
## [1] 3000
## Current iteration :
## [1] 3500
## Current iteration :
## [1] 4000
## Current iteration :
## [1] 4500
## Current iteration :
## [1] 5000
summary(model_bqr2)
##
## Type of dependent variable: continuous
## Lasso variable selection: no
## Normal approximation of posterior: yes
## Estimated quantile: 0.25
## Lower credible bound: 0.025
## Upper credible bound: 0.975
## Number of burnin draws: 0
## Number of retained draws: 5000
##
##
## Summary of the estimated beta:
##
## Bayes Estimate lower upper adj.lower adj.upper
## (Intercept) 0.251 -18.77 19.6 -39.4 39.9
## X1 158.876 152.23 164.8 -135.1 452.9
## X2 -42.729 -49.72 -32.0 -469.9 384.5
## X3 12.039 5.06 18.7 -463.1 487.1
## X4 28.289 20.16 34.7 -324.4 380.9
## X5 -17.008 -18.75 -15.2 -131.8 97.8
## X6 15.850 14.89 16.7 -12.5 44.2
#Model BayesQR (tau=0,5)
model_bqr3<-bayesQR(formula=Y~X1+X2+X3+X4+X5+X6, data=dataku, quantile=0.5, ndraw=5000)
## Current iteration :
## [1] 500
## Current iteration :
## [1] 1000
## Current iteration :
## [1] 1500
## Current iteration :
## [1] 2000
## Current iteration :
## [1] 2500
## Current iteration :
## [1] 3000
## Current iteration :
## [1] 3500
## Current iteration :
## [1] 4000
## Current iteration :
## [1] 4500
## Current iteration :
## [1] 5000
summary(model_bqr3)
##
## Type of dependent variable: continuous
## Lasso variable selection: no
## Normal approximation of posterior: yes
## Estimated quantile: 0.5
## Lower credible bound: 0.025
## Upper credible bound: 0.975
## Number of burnin draws: 0
## Number of retained draws: 5000
##
##
## Summary of the estimated beta:
##
## Bayes Estimate lower upper adj.lower adj.upper
## (Intercept) -6.92 -26.7 13.0 -43.11 29.26
## X1 270.47 266.6 272.9 123.28 417.66
## X2 -66.66 -73.0 -55.7 -376.85 243.52
## X3 22.28 20.9 23.4 -1.06 45.63
## X4 77.54 70.1 82.5 -134.53 289.61
## X5 -37.66 -38.4 -36.7 -65.35 -9.97
## X6 32.63 31.8 33.3 7.84 57.42
#Model BayesQR (tau=0,75)
model_bqr4<-bayesQR(formula=Y~X1+X2+X3+X4+X5+X6, data=dataku, quantile=0.75, ndraw=5000)
## Current iteration :
## [1] 500
## Current iteration :
## [1] 1000
## Current iteration :
## [1] 1500
## Current iteration :
## [1] 2000
## Current iteration :
## [1] 2500
## Current iteration :
## [1] 3000
## Current iteration :
## [1] 3500
## Current iteration :
## [1] 4000
## Current iteration :
## [1] 4500
## Current iteration :
## [1] 5000
summary(model_bqr4)
##
## Type of dependent variable: continuous
## Lasso variable selection: no
## Normal approximation of posterior: yes
## Estimated quantile: 0.75
## Lower credible bound: 0.025
## Upper credible bound: 0.975
## Number of burnin draws: 0
## Number of retained draws: 5000
##
##
## Summary of the estimated beta:
##
## Bayes Estimate lower upper adj.lower adj.upper
## (Intercept) -4.96 -24.1 15.0 -50.3 40.4
## X1 302.19 296.4 306.7 30.9 573.5
## X2 -32.34 -37.6 -24.5 -235.2 170.5
## X3 15.24 10.9 19.4 -213.8 244.2
## X4 115.81 102.3 124.7 -240.8 472.4
## X5 -35.19 -39.2 -32.2 -265.9 195.5
## X6 28.63 25.9 32.1 -174.5 231.8
#Model BayesQR (tau=0,9)
model_bqr5<-bayesQR(formula=Y~X1+X2+X3+X4+X5+X6, data=dataku, quantile=0.9, ndraw=5000)
## Current iteration :
## [1] 500
## Current iteration :
## [1] 1000
## Current iteration :
## [1] 1500
## Current iteration :
## [1] 2000
## Current iteration :
## [1] 2500
## Current iteration :
## [1] 3000
## Current iteration :
## [1] 3500
## Current iteration :
## [1] 4000
## Current iteration :
## [1] 4500
## Current iteration :
## [1] 5000
summary(model_bqr5)
##
## Type of dependent variable: continuous
## Lasso variable selection: no
## Normal approximation of posterior: yes
## Estimated quantile: 0.9
## Lower credible bound: 0.025
## Upper credible bound: 0.975
## Number of burnin draws: 0
## Number of retained draws: 5000
##
##
## Summary of the estimated beta:
##
## Bayes Estimate lower upper adj.lower adj.upper
## (Intercept) -0.909 -20.87 18.359 -37.1 35.2
## X1 274.022 258.40 286.883 -702.1 1250.2
## X2 -39.241 -57.26 -20.849 -699.5 621.0
## X3 -4.735 -8.73 -0.702 -203.6 194.1
## X4 83.863 68.80 100.834 -488.8 656.5
## X5 -43.699 -46.79 -40.299 -199.9 112.5
## X6 46.349 43.11 49.640 -85.4 178.1
#Menghitung MSE BQR
coef_mean_1 <- colMeans(model_bqr1[[1]]$betadraw)
coef_mean_2 <- colMeans(model_bqr2[[1]]$betadraw)
coef_mean_3 <- colMeans(model_bqr3[[1]]$betadraw)
coef_mean_4 <- colMeans(model_bqr4[[1]]$betadraw)
coef_mean_5 <- colMeans(model_bqr5[[1]]$betadraw)
coef_mean_1
## [1] -2.886398 115.806747 -26.410624 -4.572015 -9.593254 -8.434991 11.418005
coef_mean_2
## [1] 0.2505068 158.8763354 -42.7289515 12.0390369 28.2889925 -17.0082026
## [7] 15.8503161
coef_mean_3
## [1] -6.92380 270.46780 -66.66062 22.28209 77.53815 -37.65909 32.62962
coef_mean_4
## [1] -4.957659 302.194766 -32.337780 15.236409 115.814434 -35.190211 28.629332
coef_mean_5
## [1] -0.9093332 274.0215154 -39.2414172 -4.7353034 83.8629574 -43.6990743
## [7] 46.3493823
X_mat <- model.matrix(~ X1 + X2 + X3 + X4 + X5 + X6, data=dataku)
pred_1 <- X_mat %*% coef_mean_1
pred_2 <- X_mat %*% coef_mean_2
pred_3 <- X_mat %*% coef_mean_3
pred_4 <- X_mat %*% coef_mean_4
pred_5 <- X_mat %*% coef_mean_5
y <- dataku$Y
mse_1 <- mean((y - pred_1)^2)
mse_2 <- mean((y - pred_2)^2)
mse_3 <- mean((y - pred_3)^2)
mse_4 <- mean((y - pred_4)^2)
mse_5 <- mean((y - pred_5)^2)
#Print MSE BQR
mse_1
## [1] 1007627
mse_2
## [1] 860995.3
mse_3
## [1] 538859.9
mse_4
## [1] 494972.4
mse_5
## [1] 608475.7
library(quantreg)
## Loading required package: SparseM
#model_reglin
model_reglin<-lm(Y~X1+X2+X3+X4+X5+X6, data=dataku)
summary(model_reglin)
##
## Call:
## lm(formula = Y ~ X1 + X2 + X3 + X4 + X5 + X6, data = dataku)
##
## Residuals:
## Min 1Q Median 3Q Max
## -728.79 -390.74 -51.25 249.29 2169.58
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -10573.91 5020.60 -2.106 0.043391 *
## X1 537.08 69.02 7.781 8.85e-09 ***
## X2 -236.04 129.06 -1.829 0.077039 .
## X3 47.54 46.83 1.015 0.317875
## X4 -15.02 146.53 -0.103 0.918988
## X5 -20.26 37.82 -0.536 0.596049
## X6 159.60 41.98 3.802 0.000631 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 612.5 on 31 degrees of freedom
## Multiple R-squared: 0.7129, Adjusted R-squared: 0.6573
## F-statistic: 12.83 on 6 and 31 DF, p-value: 3.07e-07
#Regresi Kuantil Multivariate
model_kuantil=rq(Y~X1+X2+X3+X4+X5+X6,data=dataku,tau=c(0.1,0.25,0.5,0.75,0.9))
summary(model_kuantil,se='iid',method='br')
## Warning in rq.fit.br(x, y, tau = tau, ...): Solution may be nonunique
## Warning in rq.fit.br(x, y, tau = tau, ...): Solution may be nonunique
##
## Call: rq(formula = Y ~ X1 + X2 + X3 + X4 + X5 + X6, tau = c(0.1, 0.25,
## 0.5, 0.75, 0.9), data = dataku)
##
## tau: [1] 0.1
##
## Coefficients:
## Value Std. Error t value Pr(>|t|)
## (Intercept) -4112.96822 1904.37090 -2.15975 0.03864
## X1 161.89059 26.18182 6.18332 0.00000
## X2 37.30031 48.95222 0.76197 0.45183
## X3 28.85458 17.76459 1.62427 0.11444
## X4 7.07020 55.57879 0.12721 0.89960
## X5 6.25019 14.34451 0.43572 0.66606
## X6 41.45913 15.92209 2.60387 0.01402
##
## Call: rq(formula = Y ~ X1 + X2 + X3 + X4 + X5 + X6, tau = c(0.1, 0.25,
## 0.5, 0.75, 0.9), data = dataku)
##
## tau: [1] 0.25
##
## Coefficients:
## Value Std. Error t value Pr(>|t|)
## (Intercept) -2214.92711 1510.82767 -1.46604 0.15271
## X1 341.90618 20.77128 16.46053 0.00000
## X2 -146.84592 38.83612 -3.78117 0.00067
## X3 13.95563 14.09349 0.99022 0.32973
## X4 63.79598 44.09329 1.44684 0.15798
## X5 -18.59135 11.38018 -1.63366 0.11245
## X6 45.02571 12.63175 3.56449 0.00120
##
## Call: rq(formula = Y ~ X1 + X2 + X3 + X4 + X5 + X6, tau = c(0.1, 0.25,
## 0.5, 0.75, 0.9), data = dataku)
##
## tau: [1] 0.5
##
## Coefficients:
## Value Std. Error t value Pr(>|t|)
## (Intercept) -2892.17762 2818.17233 -1.02626 0.31271
## X1 296.20846 38.74501 7.64507 0.00000
## X2 -101.39057 72.44167 -1.39962 0.17156
## X3 25.61717 26.28883 0.97445 0.33738
## X4 106.28656 82.24796 1.29227 0.20581
## X5 -23.78510 21.22764 -1.12048 0.27112
## X6 55.98828 23.56221 2.37619 0.02385
##
## Call: rq(formula = Y ~ X1 + X2 + X3 + X4 + X5 + X6, tau = c(0.1, 0.25,
## 0.5, 0.75, 0.9), data = dataku)
##
## tau: [1] 0.75
##
## Coefficients:
## Value Std. Error t value Pr(>|t|)
## (Intercept) -2408.26875 3278.08668 -0.73466 0.46807
## X1 294.89221 45.06804 6.54327 0.00000
## X2 -14.55075 84.26385 -0.17268 0.86402
## X3 37.33494 30.57907 1.22093 0.23132
## X4 123.54339 95.67050 1.29134 0.20613
## X5 -27.78654 24.69190 -1.12533 0.26909
## X6 49.00715 27.40747 1.78809 0.08354
##
## Call: rq(formula = Y ~ X1 + X2 + X3 + X4 + X5 + X6, tau = c(0.1, 0.25,
## 0.5, 0.75, 0.9), data = dataku)
##
## tau: [1] 0.9
##
## Coefficients:
## Value Std. Error t value Pr(>|t|)
## (Intercept) -8993.86681 6415.56514 -1.40188 0.17088
## X1 807.69273 88.20297 9.15721 0.00000
## X2 -67.54223 164.91334 -0.40956 0.68494
## X3 56.13082 59.84649 0.93791 0.35554
## X4 236.41288 187.23736 1.26264 0.21614
## X5 -31.32182 48.32469 -0.64815 0.52166
## X6 126.16113 53.63935 2.35203 0.02520
#Model QR (tau=0,1)
model_qr1<-rq(formula=Y~X1+X2+X3+X4+X5+X6, data=dataku, tau=0.1)
summary(model_qr1)
##
## Call: rq(formula = Y ~ X1 + X2 + X3 + X4 + X5 + X6, tau = 0.1, data = dataku)
##
## tau: [1] 0.1
##
## Coefficients:
## coefficients lower bd upper bd
## (Intercept) -4112.96822 -21860.44668 1479.94516
## X1 161.89059 -36.82273 602.57327
## X2 37.30031 -312.82233 128.08107
## X3 28.85458 -59.53005 91.01534
## X4 7.07020 -94.08661 1151.39743
## X5 6.25019 -86.45939 51.25943
## X6 41.45913 -6.88510 223.10398
#Model QR (tau=0,25)
model_qr2<-rq(formula=Y~X1+X2+X3+X4+X5+X6, data=dataku, tau=0.25)
summary(model_qr2)
##
## Call: rq(formula = Y ~ X1 + X2 + X3 + X4 + X5 + X6, tau = 0.25, data = dataku)
##
## tau: [1] 0.25
##
## Coefficients:
## coefficients lower bd upper bd
## (Intercept) -2214.92711 -13519.92113 4675.26940
## X1 341.90618 76.86746 481.77411
## X2 -146.84592 -272.92471 110.89820
## X3 13.95563 -58.74046 81.29309
## X4 63.79598 -271.30897 452.34809
## X5 -18.59135 -65.16493 30.22540
## X6 45.02571 0.15023 134.28945
#Model QR (tau=0,5)
model_qr3<-rq(formula=Y~X1+X2+X3+X4+X5+X6, data=dataku, tau=0.5)
summary(model_qr3)
##
## Call: rq(formula = Y ~ X1 + X2 + X3 + X4 + X5 + X6, tau = 0.5, data = dataku)
##
## tau: [1] 0.5
##
## Coefficients:
## coefficients lower bd upper bd
## (Intercept) -2892.17762 -13123.98816 4901.47399
## X1 296.20846 177.94514 384.80625
## X2 -101.39057 -209.73262 -19.48271
## X3 25.61717 8.57284 54.35415
## X4 106.28656 -351.01756 161.99384
## X5 -23.78510 -54.38691 -21.93801
## X6 55.98828 17.62289 110.77815
#Model QR (tau=0,75)
model_qr4<-rq(formula=Y~X1+X2+X3+X4+X5+X6, data=dataku, tau=0.75)
summary(model_qr4)
##
## Call: rq(formula = Y ~ X1 + X2 + X3 + X4 + X5 + X6, tau = 0.75, data = dataku)
##
## tau: [1] 0.75
##
## Coefficients:
## coefficients lower bd upper bd
## (Intercept) -2408.26875 -9030.73126 5995.86363
## X1 294.89221 285.24958 728.66350
## X2 -14.55075 -281.75007 48.88820
## X3 37.33494 -33.64958 103.97997
## X4 123.54339 -454.16546 196.60930
## X5 -27.78654 -80.21199 -2.38181
## X6 49.00715 -10.24541 143.86855
#Model QR (tau=0,9)
model_qr5<-rq(formula=Y~X1+X2+X3+X4+X5+X6, data=dataku, tau=0.9)
summary(model_qr5)
##
## Call: rq(formula = Y ~ X1 + X2 + X3 + X4 + X5 + X6, tau = 0.9, data = dataku)
##
## tau: [1] 0.9
##
## Coefficients:
## coefficients lower bd upper bd
## (Intercept) -8993.86681 -17924.15389 14037.15020
## X1 807.69273 220.97357 838.94239
## X2 -67.54223 -973.25245 422.93069
## X3 56.13082 -236.32525 307.78996
## X4 236.41288 -3498.32501 422.44846
## X5 -31.32182 -106.24606 141.04407
## X6 126.16113 8.69069 374.44596
#Menghitung MSE
nilai_Y<-dataku$Y
pred_1<-predict(model_qr1)
pred_2<-predict(model_qr2)
pred_3<-predict(model_qr3)
pred_4<-predict(model_qr4)
pred_5<-predict(model_qr5)
#hitung residual
res1<-nilai_Y - pred_1
res2<-nilai_Y - pred_2
res3<-nilai_Y - pred_3
res4<-nilai_Y - pred_4
res5<-nilai_Y - pred_5
#hitung nilai MSE
n<-length(nilai_Y)
mse_qr1<-sum(res1^2)/n
mse_qr2<-sum(res2^2)/n
mse_qr3<-sum(res3^2)/n
mse_qr4<-sum(res4^2)/n
mse_qr5<-sum(res5^2)/n
#tampilkan nilai mse
mse_qr1
## [1] 887135.4
mse_qr2
## [1] 542102.6
mse_qr3
## [1] 485312.9
mse_qr4
## [1] 482436.3
mse_qr5
## [1] 1428778
# Mengambil CI 95% untuk setiap tau
CI_list <- lapply(summary(model_kuantil,se='iid',method='br'), function(m) {
coef_mat <- m$coefficients
est <- coef_mat[, 1]
se <- coef_mat[, 2]
lower <- est - 1.96 * se
upper <- est + 1.96 * se
data.frame(
Estimate = est,
CI_lower = lower,
CI_upper = upper
)
})
## Warning in rq.fit.br(x, y, tau = tau, ...): Solution may be nonunique
## Warning in rq.fit.br(x, y, tau = tau, ...): Solution may be nonunique
#Print CI 95% - QR
CI_list
## [[1]]
## Estimate CI_lower CI_upper
## (Intercept) -4112.968225 -7845.535185 -380.40126
## X1 161.890589 110.574224 213.20695
## X2 37.300310 -58.646047 133.24667
## X3 28.854582 -5.964024 63.67319
## X4 7.070200 -101.864238 116.00464
## X5 6.250189 -21.865047 34.36543
## X6 41.459126 10.251828 72.66642
##
## [[2]]
## Estimate CI_lower CI_upper
## (Intercept) -2214.92711 -5176.14935 746.295131
## X1 341.90618 301.19448 382.617887
## X2 -146.84592 -222.96471 -70.727125
## X3 13.95563 -13.66762 41.578879
## X4 63.79598 -22.62687 150.218831
## X5 -18.59135 -40.89650 3.713798
## X6 45.02571 20.26748 69.783933
##
## [[3]]
## Estimate CI_lower CI_upper
## (Intercept) -2892.17762 -8415.79539 2631.44015
## X1 296.20846 220.26823 372.14869
## X2 -101.39057 -243.37623 40.59510
## X3 25.61717 -25.90894 77.14328
## X4 106.28656 -54.91944 267.49256
## X5 -23.78510 -65.39127 17.82107
## X6 55.98828 9.80634 102.17022
##
## [[4]]
## Estimate CI_lower CI_upper
## (Intercept) -2408.26875 -8833.31863 4016.78114
## X1 294.89221 206.55884 383.22558
## X2 -14.55075 -179.70790 150.60641
## X3 37.33494 -22.60003 97.26991
## X4 123.54339 -63.97078 311.05756
## X5 -27.78654 -76.18267 20.60959
## X6 49.00715 -4.71150 102.72579
##
## [[5]]
## Estimate CI_lower CI_upper
## (Intercept) -8993.86681 -21568.37449 3580.64086
## X1 807.69273 634.81492 980.57055
## X2 -67.54223 -390.77238 255.68793
## X3 56.13082 -61.16830 173.42995
## X4 236.41288 -130.57234 603.39811
## X5 -31.32182 -126.03821 63.39457
## X6 126.16113 21.02801 231.29424
R Markdown