library(rjags)
## Loading required package: coda
## Linked to JAGS 4.3.1
## Loaded modules: basemod,bugs
library(coda)
library(ggplot2)
Load Data
data_csv <- read.csv("C:/SKRIPSI/dataset_skripsi.csv")
head(data_csv)
## kabupaten_kota TPT rata_rata_lama_sekolah IPM kepadatan_penduduk
## 1 Nias 2.31 6.14 65.15 172.46
## 2 Mandailing Natal 7.45 8.84 72.65 77.19
## 3 Tapanuli Selatan 3.49 9.51 74.58 75.34
## 4 Tapanuli Tengah 7.81 8.87 72.77 171.13
## 5 Tapanuli Utara 1.03 10.09 76.86 83.94
## 6 Toba Samosir 1.30 10.59 77.83 94.57
## Penduduk_laki_laki PDRB Angkatan_kerja persentase_penduduk_miskin
## 1 75063 4857.80 87.981 15.10
## 2 247811 18322.32 219.010 8.86
## 3 157340 18875.45 167.625 7.01
## 4 195086 12631.77 211.032 11.50
## 5 161782 10489.70 187.311 8.54
## 6 106907 9656.27 123.572 8.04
## pengeluaran_perkapita TPAK Harapan_Lama_Sekolah laju_pertumbuhan_penduduk
## 1 7.301 81.68 13.30 1.59
## 2 10.251 63.07 13.86 1.79
## 3 11.829 75.64 13.58 1.35
## 4 10.690 74.52 13.49 2.11
## 5 12.115 81.20 13.73 1.19
## 6 12.676 80.87 13.59 1.34
Statistik Deskriptif variabel Respon
c(summary(data_csv$TPT), SD = sd(data_csv$TPT))
## Min. 1st Qu. Median Mean 3rd Qu. Max. SD
## 0.450000 2.570000 4.840000 4.553030 6.240000 8.670000 2.551358
Uji Distribusi Normal
# Uji Kolmogorov-Smirnov untuk normalitas
ks_result <- ks.test(data_csv$TPT, "pnorm", mean = mean(data_csv$TPT), sd = sd(data_csv$TPT))
## Warning in ks.test.default(data_csv$TPT, "pnorm", mean = mean(data_csv$TPT), :
## ties should not be present for the one-sample Kolmogorov-Smirnov test
# Tampilkan hasil uji K-S
print(ks_result)
##
## Asymptotic one-sample Kolmogorov-Smirnov test
##
## data: data_csv$TPT
## D = 0.11097, p-value = 0.8111
## alternative hypothesis: two-sided
# Hitung mean dan standar deviasi
mean_tpt <- mean(data_csv$TPT)
sd_tpt <- sd(data_csv$TPT)
# Plot histogram
hist(data_csv$TPT, breaks = 10, probability = TRUE,
main = "Distribusi TPT dan Kurva Normal",
xlab = "TPT (%)", col = "#AED6F1", border = "white")
# Tambah kurva normal
curve(dnorm(x, mean = mean_tpt, sd = sd_tpt),
col = "#E74C3C", lwd = 2, add = TRUE)
# Garis vertikal mean dan ±1 SD
abline(v = mean_tpt, col = "#2E86C1", lwd = 2, lty = 2)
abline(v = c(mean_tpt - sd_tpt, mean_tpt + sd_tpt),
col = "#27AE60", lwd = 2, lty = 3)
# Tentukan posisi teks pojok kanan atas
text_x <- max(data_csv$TPT)*0.95
text_y <- max(density(data_csv$TPT)$y) * 0.95
# Tambah informasi Mean & SD
text(text_x, text_y,
labels = paste("\n",
"μ =", round(mean_tpt, 2),
"\nσ =", round(sd_tpt, 2)),
col = "black", cex = 0.9, adj = c(1, 1))
# Legenda
legend("topright", legend = c("Kurva Normal", "μ (Mean)", "±1SD"),
col = c("#E74C3C", "#2E86C1", "#27AE60"), lty = c(1, 2, 3), bty = "n", cex = 0.9)
Korelasi
# Daftar nama variabel prediktor
var_names <- c(
"rata_rata_lama_sekolah", "IPM", "kepadatan_penduduk",
"PDRB", "persentase_penduduk_miskin",
"pengeluaran_perkapita", "TPAK",
"Harapan_Lama_Sekolah", "laju_pertumbuhan_penduduk"
)
# Korelasi langsung dari data asli
results_raw <- data.frame()
for (v in var_names) {
test <- cor.test(data_csv$TPT, data_csv[[v]])
results_raw <- rbind(results_raw, data.frame(
Variable = v,
Correlation = round(test$estimate, 3),
P_Value = round(test$p.value, 4),
Significant = ifelse(test$p.value < 0.05, "Yes", "No")
))
}
print(results_raw)
## Variable Correlation P_Value Significant
## cor rata_rata_lama_sekolah 0.370 0.0340 Yes
## cor1 IPM 0.478 0.0049 Yes
## cor2 kepadatan_penduduk 0.484 0.0043 Yes
## cor3 PDRB 0.468 0.0060 Yes
## cor4 persentase_penduduk_miskin -0.369 0.0346 Yes
## cor5 pengeluaran_perkapita 0.550 0.0009 Yes
## cor6 TPAK -0.794 0.0000 Yes
## cor7 Harapan_Lama_Sekolah 0.304 0.0856 No
## cor8 laju_pertumbuhan_penduduk -0.272 0.1260 No
Standarisasi Prediktor
# Standarisasi variabel
scaled_vars <- as.data.frame(scale(data_csv[, var_names[1:7]]))
colnames(scaled_vars) <- paste0("x", 1:7)
# Gabungkan dengan TPT (y)
data <- cbind(y = data_csv$TPT, scaled_vars)
# Tampilkan hanya y dan x1–x7
print(data[, c("y", paste0("x", 1:7))])
## y x1 x2 x3 x4 x5 x6
## 1 2.31 -2.41949124 -2.10515978 -0.4444385 -0.48488983 1.2344080 -1.71864012
## 2 7.45 -0.40738867 -0.32362711 -0.4873176 -0.24165819 -0.2541683 -0.31419942
## 3 3.49 0.09191086 0.13482063 -0.4881502 -0.23166610 -0.6954931 0.43705733
## 4 7.81 -0.38503197 -0.29512258 -0.4450371 -0.34445590 0.3756140 -0.10519960
## 5 1.03 0.52414030 0.67640657 -0.4842795 -0.38315161 -0.3305056 0.57321666
## 6 1.30 0.89675189 0.90681812 -0.4794952 -0.39820722 -0.4497826 0.84029843
## 7 5.99 0.07700639 0.16332516 -0.4375523 0.25404044 -0.4617102 0.36136035
## 8 6.12 -0.41484090 -0.10034168 -0.4240769 0.35374674 -0.4092284 0.42087055
## 9 5.35 0.24840773 0.27496787 -0.4191755 0.37103604 -0.4903367 0.39754256
## 10 1.23 0.36764343 0.27734325 -0.4506496 -0.36217044 -0.5857583 0.02762716
## 11 2.63 0.47942691 0.68115732 -0.4350813 -0.09557582 -0.4640958 0.88933484
## 12 8.62 0.66573270 0.74766787 -0.1648715 1.94010510 -1.5471305 0.94217989
## 13 6.33 -0.48936322 0.07543621 -0.4429892 0.44243019 -0.1659034 0.34326925
## 14 3.48 -2.16611536 -2.14554118 -0.4527334 -0.41774033 1.5421426 -1.71959228
## 15 0.84 0.46452245 -0.36163314 -0.4824252 -0.43435853 -0.2947225 -1.15924425
## 16 0.45 0.16643318 -0.40676530 -0.5035339 -0.54446753 -0.5690595 -1.02213275
## 17 1.03 0.06210193 -0.25711655 -0.4872545 -0.47047515 0.4137826 -0.83455660
## 18 4.97 -0.39993644 -0.14547384 -0.3623530 0.10539722 -0.5929149 0.37326239
## 19 5.88 -0.66076455 -0.34500550 -0.3004220 0.21705129 0.3469875 0.01048822
## 20 4.42 0.12171979 -0.13359695 -0.4900495 -0.28681888 -0.2708671 -0.14090572
## 21 5.75 0.03229300 -0.44002057 -0.4897750 -0.28131440 -0.4855656 -0.72172526
## 22 3.43 -0.34031858 0.05168244 -0.4728700 0.10630262 -0.4450115 0.49466320
## 23 4.84 -0.38503197 0.34147842 -0.4721499 0.02838244 -0.2016865 0.72270628
## 24 2.57 -1.89038279 -2.03627385 -0.4643725 -0.49175312 2.8303337 -1.96286997
## 25 0.80 -1.72643369 -2.21680249 -0.4277720 -0.53043601 3.0736587 -2.15615910
## 26 6.79 0.78496841 0.72628948 3.0777540 -0.44942895 0.3565297 0.65415053
## 27 4.47 0.21859880 0.33435229 0.8688713 -0.37014887 0.5449872 0.40087513
## 28 8.62 1.63452283 1.53154225 1.1204343 -0.27949042 -0.6406257 0.98693156
## 29 6.24 1.09796214 0.98758094 1.5574436 -0.43845486 -0.1038794 1.17784028
## 30 8.67 1.66433176 2.04224828 3.4846215 4.90657407 -0.4593247 2.26759104
## 31 6.10 1.34388579 0.97332868 0.9523026 -0.31143718 -1.2250828 0.31232395
## 32 7.57 1.29172017 0.97095330 0.1562243 -0.42941062 -0.7336617 0.30518273
## 33 3.67 -0.54898107 -0.58491856 -0.2088263 -0.44755619 1.1580708 -1.08354727
## x7
## 1 1.00495590
## 2 -1.49524084
## 3 0.19350032
## 4 0.04303174
## 5 0.94046936
## 6 0.89613487
## 7 -1.25744674
## 8 -0.62735954
## 9 -0.27537054
## 10 1.42546185
## 11 1.29783225
## 12 -0.88396222
## 13 0.22440012
## 14 0.68386669
## 15 1.67534718
## 16 1.57861737
## 17 1.70490351
## 18 -0.87993181
## 19 -0.16923645
## 20 -0.10609338
## 21 0.52936769
## 22 -0.13430624
## 23 -1.27894225
## 24 0.94987365
## 25 0.92434773
## 26 -0.40568708
## 27 -0.59242934
## 28 -0.43524341
## 29 -1.14190836
## 30 -1.28028572
## 31 -1.53285799
## 32 -0.71199812
## 33 -0.86381018
Model Hierarchical Bayesian
model_file <- "model_jags.txt"
writeLines(c(
"model {",
" for (i in 1:N) {",
" y[i] ~ dnorm(theta[i], sigma_y)",
" theta[i] <- beta0 + beta1 * x1[i] + beta2 * x2[i] + beta3 * x3[i] +",
" beta4 * x4[i] + beta5 * x5[i] + beta6 * x6[i] + beta7 * x7[i] + v[i]",
" v[i] ~ dnorm(0, sigma_v)",
" }",
" beta0 ~ dnorm(0, 1)",
" beta1 ~ dnorm(0, 1)",
" beta2 ~ dnorm(0, 1)",
" beta3 ~ dnorm(0, 1)",
" beta4 ~ dnorm(0, 1)",
" beta5 ~ dnorm(0, 1)",
" beta6 ~ dnorm(0, 1)",
" beta7 ~ dnorm(0, 1)",
" tau_y ~ dgamma(2, 0.1)",
" tau_v ~ dgamma(2, 0.1)",
" sigma_y <- 1 / tau_y",
" sigma_v <- 1 / tau_v",
"}"
), con = model_file)
inits <- function(chain) {
list(
beta0 = 0, beta1 = 0, beta2 = 0, beta3 = 0,
beta4 = 0, beta5 = 0, beta6 = 0, beta7 = 0,
tau_y = 1, tau_v = 1, v = rep(0, nrow(data)),
.RNG.name = "base::Wichmann-Hill",
.RNG.seed = 1234 + chain
)
}
data_jags <- list(N = nrow(data), y = data$y,
x1 = data$x1, x2 = data$x2, x3 = data$x3, x4 = data$x4,
x5 = data$x5, x6 = data$x6, x7 = data$x7)
model <- jags.model(model_file, data = data_jags, inits = inits, n.chains = 3)
## Compiling model graph
## Resolving undeclared variables
## Allocating nodes
## Graph information:
## Observed stochastic nodes: 33
## Unobserved stochastic nodes: 43
## Total graph size: 577
##
## Initializing model
update(model, 5000)
samples <- coda.samples(model, variable.names = c("beta0", "beta1", "beta2", "beta3",
"beta4", "beta5", "beta6", "beta7",
"tau_v", "tau_y", "v"),
n.iter = 50000, thin = 25)
# Ringkasan hasil
summary(samples)
##
## Iterations = 6025:56000
## Thinning interval = 25
## Number of chains = 3
## Sample size per chain = 2000
##
## 1. Empirical mean and standard deviation for each variable,
## plus standard error of the mean:
##
## Mean SD Naive SE Time-series SE
## beta0 4.063897 0.3661 0.004726 0.004828
## beta1 -0.113621 0.6672 0.008613 0.008853
## beta2 -0.010618 0.8085 0.010437 0.010968
## beta3 0.400618 0.4753 0.006136 0.006137
## beta4 0.314567 0.4020 0.005190 0.005054
## beta5 -0.345106 0.5416 0.006993 0.006993
## beta6 0.099047 0.6530 0.008430 0.008971
## beta7 -1.440665 0.4017 0.005186 0.005186
## tau_v 2.026930 1.2460 0.016085 0.018137
## tau_y 2.034944 1.2307 0.015889 0.017055
## v[1] 0.169359 1.0415 0.013446 0.013446
## v[2] 0.694685 1.0612 0.013700 0.013893
## v[3] -0.127381 0.9393 0.012126 0.012258
## v[4] 2.083585 1.3291 0.017158 0.018743
## v[5] -0.747950 1.0293 0.013288 0.013288
## v[6] -0.631289 1.0104 0.013044 0.013491
## v[7] 0.008891 0.9647 0.012455 0.011950
## v[8] 0.492746 0.9804 0.012657 0.012661
## v[9] 0.377950 0.9497 0.012261 0.012261
## v[10] -0.311050 1.0086 0.013022 0.013023
## v[11] 0.227660 1.0176 0.013137 0.013138
## v[12] 1.117287 1.1469 0.014806 0.015297
## v[13] 1.252947 1.1066 0.014286 0.014511
## v[14] 0.577580 1.0628 0.013721 0.013722
## v[15] -0.204232 1.0505 0.013561 0.013716
## v[16] -0.516411 1.0439 0.013477 0.013473
## v[17] -0.008489 0.9923 0.012811 0.012811
## v[18] -0.275783 0.9726 0.012556 0.012468
## v[19] 0.863660 1.0328 0.013333 0.012253
## v[20] 0.202018 0.9254 0.011947 0.012017
## v[21] 1.338617 1.1390 0.014705 0.015299
## v[22] -0.449672 0.9644 0.012451 0.012160
## v[23] -0.534935 1.0278 0.013268 0.013054
## v[24] 0.583744 1.0930 0.014111 0.014479
## v[25] -0.296877 1.0907 0.014081 0.014533
## v[26] 0.601673 1.1847 0.015294 0.015982
## v[27] -0.247343 0.9743 0.012578 0.013477
## v[28] 1.710141 1.2602 0.016269 0.017127
## v[29] 0.022524 1.0394 0.013419 0.013356
## v[30] -0.176189 1.3264 0.017124 0.016392
## v[31] -0.380716 1.0687 0.013797 0.014108
## v[32] 1.213050 1.1291 0.014577 0.015114
## v[33] -0.456099 1.0148 0.013101 0.012453
##
## 2. Quantiles for each variable:
##
## 2.5% 25% 50% 75% 97.5%
## beta0 3.2748 3.84132 4.085776 4.30938 4.7226
## beta1 -1.4104 -0.56648 -0.108103 0.33818 1.2035
## beta2 -1.6006 -0.54792 -0.009455 0.54487 1.5730
## beta3 -0.5253 0.08520 0.403562 0.70988 1.3324
## beta4 -0.4805 0.05375 0.311711 0.57877 1.1091
## beta5 -1.3930 -0.71252 -0.339268 0.01039 0.7305
## beta6 -1.1363 -0.33849 0.086141 0.53515 1.4110
## beta7 -2.2222 -1.70206 -1.447441 -1.18430 -0.6155
## tau_v 0.3198 1.14084 1.812053 2.64572 5.0750
## tau_y 0.3526 1.15261 1.838655 2.62391 4.9259
## v[1] -1.8508 -0.49234 0.132584 0.81116 2.3704
## v[2] -1.3124 -0.01325 0.654362 1.37858 2.9000
## v[3] -1.9380 -0.72807 -0.129844 0.46361 1.7437
## v[4] -0.2683 1.07077 2.047177 3.03767 4.6192
## v[5] -2.8104 -1.41301 -0.718819 -0.06236 1.2443
## v[6] -2.7150 -1.27682 -0.606806 0.03164 1.3181
## v[7] -1.8526 -0.62072 -0.005681 0.61958 1.9851
## v[8] -1.3816 -0.16052 0.468723 1.11731 2.4675
## v[9] -1.4501 -0.22906 0.351104 0.95827 2.2928
## v[10] -2.3189 -0.96852 -0.280216 0.33466 1.6815
## v[11] -1.7718 -0.43120 0.224968 0.87624 2.3183
## v[12] -0.9092 0.31318 1.050714 1.84841 3.5642
## v[13] -0.7409 0.47867 1.210956 1.99273 3.4884
## v[14] -1.4653 -0.11191 0.523557 1.24386 2.7718
## v[15] -2.2568 -0.90077 -0.206681 0.45815 1.9081
## v[16] -2.6897 -1.16330 -0.493626 0.11863 1.5696
## v[17] -1.9681 -0.64304 -0.011216 0.60127 1.9955
## v[18] -2.2047 -0.89879 -0.283970 0.34809 1.6889
## v[19] -1.0862 0.16950 0.822224 1.53226 2.9618
## v[20] -1.5615 -0.39726 0.181296 0.76721 2.0885
## v[21] -0.7103 0.51926 1.296649 2.11609 3.6085
## v[22] -2.3472 -1.06165 -0.445805 0.16860 1.4580
## v[23] -2.5356 -1.20310 -0.525298 0.14478 1.5373
## v[24] -1.5046 -0.13747 0.537959 1.28238 2.8725
## v[25] -2.5192 -0.98545 -0.292630 0.39071 1.8987
## v[26] -1.5866 -0.18386 0.538224 1.31368 3.1347
## v[27] -2.1258 -0.89640 -0.246187 0.37523 1.6951
## v[28] -0.5563 0.79650 1.660105 2.55508 4.2678
## v[29] -2.0175 -0.64459 0.011430 0.67808 2.1658
## v[30] -2.8932 -0.98650 -0.150367 0.62267 2.4907
## v[31] -2.5053 -1.06168 -0.368702 0.30320 1.7854
## v[32] -0.8681 0.42215 1.166030 1.96084 3.5528
## v[33] -2.4801 -1.10709 -0.449124 0.19077 1.5552
# Trace plot dan ACF
par(mfrow = c(4, 2), mar = c(2, 2, 2, 2))
plot(samples, ask = FALSE, col = "red")
autocorr.plot(samples, col = "red")
# Ergodic Mean (Cumulative Mean)
samples_matrix <- as.matrix(samples)
# Fungsi plot cumulative mean
plot_parameter_convergence <- function(param_name, samples_matrix) {
param_samples <- samples_matrix[, param_name]
plot(cumsum(param_samples) / seq_along(param_samples),
type = "l",
col = "red",
main = paste("Ergodic Mean for", param_name),
ylab = "Cumulative Mean",
xlab = "Iteration")
}
# Plot untuk beta dan precision
params <- colnames(samples_matrix)
par(mfrow = c(3, 2))
for (param in params) {
plot_parameter_convergence(param, samples_matrix)
}
# Plot ergodic mean untuk v[i]
num_v <- sum(grepl("^v\\[", colnames(samples_matrix)))
for (i in 1:num_v) {
param_name <- paste0("v[", i, "]")
plot_parameter_convergence(param_name, samples_matrix)
}
# Ringkasan hasil MCMC singkat
stat <- summary(samples)$statistics
quant <- summary(samples)$quantiles
table_result <- data.frame(
Parameter = rownames(stat),
Mean = round(stat[, "Mean"], 4),
SD = round(stat[, "SD"], 4),
MC_Error = round(stat[, "Time-series SE"], 4),
CI_2.5 = round(quant[, "2.5%"], 4),
CI_97.5 = round(quant[, "97.5%"], 4),
ESS = round(effectiveSize(samples), 2),
Rhat = if (inherits(samples, "mcmc.list"))
round(gelman.diag(samples, autoburnin = FALSE)$psrf[, "Point est."], 3) else NA
)
print(table_result)
## Parameter Mean SD MC_Error CI_2.5 CI_97.5 ESS Rhat
## beta0 beta0 4.0639 0.3661 0.0048 3.2748 4.7226 5785.98 1.001
## beta1 beta1 -0.1136 0.6672 0.0089 -1.4104 1.2035 5682.09 1.000
## beta2 beta2 -0.0106 0.8085 0.0110 -1.6006 1.5730 5499.37 1.001
## beta3 beta3 0.4006 0.4753 0.0061 -0.5253 1.3324 6000.00 1.000
## beta4 beta4 0.3146 0.4020 0.0051 -0.4805 1.1091 6349.28 1.000
## beta5 beta5 -0.3451 0.5416 0.0070 -1.3930 0.7305 6000.00 1.000
## beta6 beta6 0.0990 0.6530 0.0090 -1.1363 1.4110 5450.44 1.000
## beta7 beta7 -1.4407 0.4017 0.0052 -2.2222 -0.6155 6000.00 1.000
## tau_v tau_v 2.0269 1.2460 0.0181 0.3198 5.0750 4727.87 1.001
## tau_y tau_y 2.0349 1.2307 0.0171 0.3526 4.9259 5212.24 1.000
## v[1] v[1] 0.1694 1.0415 0.0134 -1.8508 2.3704 6000.00 1.000
## v[2] v[2] 0.6947 1.0612 0.0139 -1.3124 2.9000 5845.53 1.000
## v[3] v[3] -0.1274 0.9393 0.0123 -1.9380 1.7437 5875.00 1.000
## v[4] v[4] 2.0836 1.3291 0.0187 -0.2683 4.6192 5038.94 1.000
## v[5] v[5] -0.7480 1.0293 0.0133 -2.8104 1.2443 6000.00 1.001
## v[6] v[6] -0.6313 1.0104 0.0135 -2.7150 1.3181 5609.58 1.001
## v[7] v[7] 0.0089 0.9647 0.0120 -1.8526 1.9851 6530.39 1.000
## v[8] v[8] 0.4927 0.9804 0.0127 -1.3816 2.4675 6012.21 1.000
## v[9] v[9] 0.3780 0.9497 0.0123 -1.4501 2.2928 6000.00 1.000
## v[10] v[10] -0.3110 1.0086 0.0130 -2.3189 1.6815 6000.00 1.000
## v[11] v[11] 0.2277 1.0176 0.0131 -1.7718 2.3183 6000.00 1.000
## v[12] v[12] 1.1173 1.1469 0.0153 -0.9092 3.5642 5629.45 1.000
## v[13] v[13] 1.2529 1.1066 0.0145 -0.7409 3.4884 5830.40 1.000
## v[14] v[14] 0.5776 1.0628 0.0137 -1.4653 2.7718 6000.00 1.000
## v[15] v[15] -0.2042 1.0505 0.0137 -2.2568 1.9081 5869.97 1.001
## v[16] v[16] -0.5164 1.0439 0.0135 -2.6897 1.5696 6000.00 1.001
## v[17] v[17] -0.0085 0.9923 0.0128 -1.9681 1.9955 6000.00 1.001
## v[18] v[18] -0.2758 0.9726 0.0125 -2.2047 1.6889 6087.88 1.000
## v[19] v[19] 0.8637 1.0328 0.0123 -1.0862 2.9618 7507.47 1.000
## v[20] v[20] 0.2020 0.9254 0.0120 -1.5615 2.0885 5935.49 1.000
## v[21] v[21] 1.3386 1.1390 0.0153 -0.7103 3.6085 5586.20 1.000
## v[22] v[22] -0.4497 0.9644 0.0122 -2.3472 1.4580 6323.87 1.000
## v[23] v[23] -0.5349 1.0278 0.0131 -2.5356 1.5373 6215.71 1.000
## v[24] v[24] 0.5837 1.0930 0.0145 -1.5046 2.8725 5718.88 1.000
## v[25] v[25] -0.2969 1.0907 0.0145 -2.5192 1.8987 5650.54 1.000
## v[26] v[26] 0.6017 1.1847 0.0160 -1.5866 3.1347 5538.58 1.000
## v[27] v[27] -0.2473 0.9743 0.0135 -2.1258 1.6951 5287.68 1.000
## v[28] v[28] 1.7101 1.2602 0.0171 -0.5563 4.2678 5420.31 1.000
## v[29] v[29] 0.0225 1.0394 0.0134 -2.0175 2.1658 6058.21 1.000
## v[30] v[30] -0.1762 1.3264 0.0164 -2.8932 2.4907 6582.27 1.000
## v[31] v[31] -0.3807 1.0687 0.0141 -2.5053 1.7854 5747.29 1.000
## v[32] v[32] 1.2131 1.1291 0.0151 -0.8681 3.5528 5604.28 1.001
## v[33] v[33] -0.4561 1.0148 0.0125 -2.4801 1.5552 6671.20 1.000
Ringkasan dan Estimasi
# Tambahkan ke dalam data
data$kabupaten <- data_csv$kabupaten_kota
# Hitung ulang theta_hat
summary_samples <- summary(samples)
beta_means <- summary_samples$statistics[, "Mean"]
v_i <- beta_means[grep("^v\\[", names(beta_means))]
beta_values <- beta_means[paste0("beta", 0:7)]
data$theta_hat <- beta_values["beta0"] +
beta_values["beta7"] * data$x7 +
v_i
# Tampilkan hasil
print(data[, c("kabupaten", "theta_hat")])
## kabupaten theta_hat
## 1 Nias 2.785451
## 2 Mandailing Natal 6.912722
## 3 Tapanuli Selatan 3.657747
## 4 Tapanuli Tengah 6.085487
## 5 Tapanuli Utara 1.961046
## 6 Toba Samosir 2.141578
## 7 Labuhan Batu 5.884347
## 8 Asahan 5.460457
## 9 Simalungun 4.838563
## 10 Dairi 1.699234
## 11 Karo 2.421815
## 12 Deli Serdang 6.454677
## 13 Langkat 4.993559
## 14 Nias Selatan 3.656254
## 15 Humbang Hasundutan 1.446051
## 16 Pakpak Bharat 1.273228
## 17 Samosir 1.599213
## 18 Serdang Bedagai 5.055800
## 19 Batu Bara 5.171370
## 20 Padang Lawas Utara 4.418759
## 21 Padang Lawas 4.639872
## 22 Labuhan Batu Selatan 3.807715
## 23 Labuhan Batu Utara 5.371489
## 24 Nias Utara 3.279192
## 25 Nias Barat 2.435345
## 26 Kota Sibolga 5.250029
## 27 Kota Tanjung Balai 4.670045
## 28 Kota Pematang Siantar 6.401078
## 29 Kota Tebing Tinggi 5.731528
## 30 Kota Medan 5.732170
## 31 Kota Binjai 5.891515
## 32 Kota Padangsidimpuan 6.302698
## 33 Kota Gunungsitoli 4.852258
Mean Squared Error (MSE)
MSE <- mean((data$y - data$theta_hat)^2)
print(paste("MSE:", round(MSE, 4)))
## [1] "MSE: 1.1683"