Tahapan Pengestimasian Direct dan Traditional Indirect Estimation:
library(readxl)
Jatim_pengeluaran <- read_excel("D:/Kuliah/tkt IV/The Journey Start/SAE/Data SAE/Jatim_pengeluaran.xlsx",
col_types = c("text", "text", "numeric",
"text", "text", "numeric", "numeric",
"numeric", "numeric"))
Xnew_ya <- read_excel("D:/Kuliah/tkt IV/The Journey Start/SAE/Data SAE/Xnew ya.xlsx",
sheet = "Sheet2", col_types = c("text", "text", "text", "text", "text", "text",
"numeric", "numeric", "numeric", "numeric", "numeric", "numeric",
"numeric", "numeric", "numeric", "numeric", "numeric", "numeric",
"numeric", "numeric", "numeric", "numeric", "numeric", "numeric",
"numeric", "numeric", "numeric", "numeric", "numeric", "numeric",
"numeric", "numeric", "numeric", "numeric", "numeric", "numeric",
"numeric", "numeric", "numeric", "numeric", "numeric"))
muatan <- read_excel("D:/Kuliah/tkt IV/The Journey Start/SAE/Data SAE/muatan.xlsx",
sheet = "Sheet3", col_types = c("text","text", "numeric"))
#variabel_X_signifikan
x1 <- aggregate(Xnew_ya$X1, list(Xnew_ya$R102),sum,na.rm=TRUE)
x12 <- aggregate(Xnew_ya$X12, list(Xnew_ya$R102),sum,na.rm=TRUE)
x13 <- aggregate(Xnew_ya$X13, list(Xnew_ya$R102),sum,na.rm=TRUE)
x20 <- aggregate(Xnew_ya$X20, list(Xnew_ya$R102),sum,na.rm=TRUE)
x27 <- aggregate(Xnew_ya$X27, list(Xnew_ya$R102),sum,na.rm=TRUE)
x31 <- aggregate(Xnew_ya$X31, list(Xnew_ya$R102),sum,na.rm=TRUE)
dt1 <- data.frame(cbind("bs_o"=Jatim_pengeluaran$`No Kode Sampel`,"pop"=Jatim_pengeluaran$Exp_krb,"kab"=Jatim_pengeluaran$R102))
mat_x <- matrix(cbind(1,x1$x,x12$x,x13$x,x20$x,x27$x,x31$x),ncol=7)
#Jumlah ruta tiap kab (Ni)
ht_rt <- function(dt1){
r102 = unique(dt1$kab)
de = data.frame(r102)
Peng_jatim <- data.frame(dt1)
ds <- list()
for (i in r102){
ds[[i]] <- subset(Peng_jatim, kab %in% i)
}
Ni <- c()
for(k in r102){
Ni[k] <- length(ds[[k]]$kab)
}
N_desa <- data.frame("desa"=Ni)
return(N_desa)
}
Ni_ruta <- ht_rt(dt1)Fungsi Estimasi untuk metode Two stage cluster SRSWOR-SRSWOR
resampling_2st_srs = function(dataY, dataX, f1, f2, iterasi ){
est_2st_cl <- function(dataY, f1,f2){
sp_2_ct <- function(dataY,f1,f2){
samp_th1 <- function(mydata,f1){ #pengambilan sampel tahap 1 secara SRSWOR
bs_u <-unique(mydata$bs_o)
bs_u <- data.frame(bs_u)
samp_bs <- sample(1:nrow(bs_u), size = f1*nrow(bs_u), replace = FALSE)
sampi <- bs_u[samp_bs,]
sampi2 <- as.data.frame(paste(sampi))
colnames(sampi2) <- c("sampel_bs")
return(sampi2)
}
hs_sp_t1 <- samp_th1(dataY,f1)
samp_th2 <- function(dt_kb_1,hs_sp_1,f2){ #Pengambilan sampel tahap 2 secara SRSWOR
index <- c(1:nrow(hs_sp_t1))
hs_sp_1 <- cbind(index,hs_sp_t1)
dt_2 <- list();dt_3 <- list() ; samp_rt <- list()
hs_samp_srs <- list() ; hs_smpl_srs <- list() ; Mi <- list()
for(i in hs_sp_1$index){ #Pengambilan sampel ssu scr srs wor
dt_2[[i]] <- subset(hs_sp_1, index %in% i)
dt_3[[i]] <- subset(dt_kb_1, bs_o %in% dt_2[[i]]$sampel_bs)
Mi[[i]] <- nrow(dt_3[[i]])
samp_rt[[i]] <- sample(1:length(dt_3[[i]]$pop), size = f2*length(dt_3[[i]]$pop), replace = FALSE)
hs_samp_srs[[i]] <- dt_3[[i]][samp_rt[[i]],]
hs_smpl_srs[[i]] <- cbind(hs_samp_srs[[i]],Mi[[i]])
}
gb_samp = data.frame()
for (j in hs_sp_1$index) {
gb_samp <- rbind(gb_samp,hs_smpl_srs[[j]])
}
return(gb_samp)
}
#Estimasi DIRECT
res_sp <- samp_th2(dataY,hs_sp_t1,f2)
res_sp$pop <- as.numeric(paste(res_sp$pop))
res_sp$bs_o <- as.character(paste(res_sp$bs_o))
res_sp$`Mi[[i]]` <- as.numeric(paste(res_sp$`Mi[[i]]`))
dataY$bs_o <- as.character(paste(dataY$bs_o))
ybar_i. <- aggregate(res_sp$pop, list(res_sp$bs_o), mean)
Mi <- aggregate(res_sp$`Mi[[i]]`, list(res_sp$bs_o), mean)
mi <- Mi$x*f2
n <- length(unique(res_sp$bs_o))
N <- length(unique(dataY$bs_o))
Y_cap_est <- (N/n) * (sum(Mi$x * ybar_i.$x))
ay <- rep(ybar_i.$x,length(res_sp$pop))
sw_i_2 <- 1/(mi-1) * sum((res_sp$pop - ay)^2)
sb_2 <- 1/(n-1) * sum((Mi$x*ybar_i.$x) - (1/n *(sum(Mi$x*ybar_i.$x)))^2)
var_cap_est <- (N^2 * (1/n - 1/N) * sb_2 ) + (N/n * sum(Mi$x^2 * (1/mi - 1/Mi$x) * sw_i_2))
rel_var_y <- var_cap_est/Y_cap_est
hasil <- cbind(Y_cap_est,var_cap_est,rel_var_y)
return(hasil)
}
r102 = unique(dt2$kab)
datakoeh = list()
for(b in r102){
datakoeh[[b]] <- subset(dt2, kab %in% b)
}
try_ye <- list()
for(c in r102){
try_ye[[c]] <- sp_2_ct(data = datakoeh[[c]],f1,f2)
}
this_it = data.frame()
for (q in r102) {
this_it <- rbind(this_it, try_ye[[q]])
}
this_it = cbind(r102, this_it)
}
#fungsi estimasi synthetic dan composite
Y_synth_comp <- function(x,y,vary){
#synthetic function
xty <- t(x)%*%y
xtx <- t(x)%*%x
bt_xy <- solve(xtx)%*%xty
bt_xy = t(bt_xy)
Y_syn <- c(rep(0,nrow(x)))
for(i in 1: ncol(x)){
Y_syn = Y_syn + x[,i]*bt_xy[,i]
}
Y_syn <- as.data.frame(Y_syn)
mse_syn <- ((Y_syn-y)^2)-vary
rel_mse_syn <- mse_syn/Y_syn
#composite function
tetha_i <- mse_syn/(vary+mse_syn)
tetha_i[tetha_i<0]<-0
tetha_i[tetha_i>1]<-1
Y_comp <- (tetha_i*y) + ((1-tetha_i)*Y_syn)
mse_comp <- (1-tetha_i)*mse_syn
rel_mse_comp <- mse_comp/Y_comp
Syn_Comp <- data.frame(cbind(Y_syn, Y_comp,mse_syn, mse_comp, rel_mse_syn, rel_mse_comp))
colnames(Syn_Comp) = c("Y syn","Y comp", "MSE syn", "MSE comp", "Relatif MSE syn","Relatif MSE comp")
return(Syn_Comp)
}
dataset = list()
est_dataset = list()
for( j in 1:iterasi){
set.seed(j)
dataset[[j]] = est_2st_cl(dataY, f1, f2)
est_dataset [[j]] = Y_synth_comp(dataX,dataset[[j]]$Y_cap_est,dataset[[j]]$var_cap_est)
}
databind = data.frame()
databind_ind = data.frame()
for (b in 1:iterasi) {
databind = rbind(databind, dataset[[b]])
databind_ind = rbind(databind_ind, est_dataset[[b]])
}
datafinal = aggregate(list(databind$Y_cap_est,databind_ind$`Y syn`,databind_ind$`Y comp`,
databind$var_cap_est,databind_ind$`MSE syn`,databind_ind$`MSE comp`
,databind$rel_var_y,databind_ind$`Relatif MSE syn`,databind_ind$`Relatif MSE comp`),
by=list(databind$r102),
FUN=mean)
colnames(datafinal) = c("kab", "Y dir", "Y syn","Y com", "MSE dir","MSE syn","MSE com",
"Relatif MSE dir","Relatif MSE syn","Relatif MSE com")
datafinal$Efisiensi_syn <- datafinal$`MSE syn`/datafinal$`MSE dir`
datafinal$Efisiensi_com <- datafinal$`MSE com`/datafinal$`MSE dir`
return(datafinal)
}b. MEnggunakan rumus MSE synthetic yang perbaikanresampling_2st_srs_2 = function(dataY, dataX, Ni, f1, f2, iterasi ){
est_2st_cl <- function(dataY, f1,f2){
sp_2_ct <- function(dataY,f1,f2){
samp_th1 <- function(mydata,f1){ #pengambilan sampel tahap 1
bs_u <-unique(mydata$bs_o)
bs_u <- data.frame(bs_u)
samp_bs <- sample(1:nrow(bs_u), size = f1*nrow(bs_u), replace = FALSE)
sampi <- bs_u[samp_bs,]
sampi2 <- as.data.frame(paste(sampi))
colnames(sampi2) <- c("sampel_bs")
return(sampi2)
}
hs_sp_t1 <- samp_th1(dataY,f1)
samp_th2 <- function(dt_kb_1,hs_sp_1,f2){
index <- c(1:nrow(hs_sp_t1))
hs_sp_1 <- cbind(index,hs_sp_t1)
dt_2 <- list();dt_3 <- list() ; samp_rt <- list()
hs_samp_srs <- list() ; hs_smpl_srs <- list() ; Mi <- list()
for(i in hs_sp_1$index){ #Pengambilan sampel ssu scr srs wor
dt_2[[i]] <- subset(hs_sp_1, index %in% i)
dt_3[[i]] <- subset(dt_kb_1, bs_o %in% dt_2[[i]]$sampel_bs)
Mi[[i]] <- nrow(dt_3[[i]])
samp_rt[[i]] <- sample(1:length(dt_3[[i]]$pop), size = f2*length(dt_3[[i]]$pop), replace = FALSE)
hs_samp_srs[[i]] <- dt_3[[i]][samp_rt[[i]],]
hs_smpl_srs[[i]] <- cbind(hs_samp_srs[[i]],Mi[[i]])
}
gb_samp = data.frame()
for (j in hs_sp_1$index) {
gb_samp <- rbind(gb_samp,hs_smpl_srs[[j]])
}
return(gb_samp)
}
res_sp <- samp_th2(dataY,hs_sp_t1,f2)
res_sp$pop <- as.numeric(paste(res_sp$pop))
res_sp$bs_o <- as.character(paste(res_sp$bs_o))
res_sp$`Mi[[i]]` <- as.numeric(paste(res_sp$`Mi[[i]]`))
dataY$bs_o <- as.character(paste(dataY$bs_o))
ybar_i. <- aggregate(res_sp$pop, list(res_sp$bs_o), mean)
Mi <- aggregate(res_sp$`Mi[[i]]`, list(res_sp$bs_o), mean)
mi <- Mi$x*f2
n <- length(unique(res_sp$bs_o))
N <- length(unique(dataY$bs_o))
Y_cap_est <- (N/n) * (sum(Mi$x * ybar_i.$x))
ay <- rep(ybar_i.$x,length(res_sp$pop))
sw_i_2 <- 1/(mi-1) * sum((res_sp$pop - ay)^2)
sb_2 <- 1/(n-1) * sum((Mi$x*ybar_i.$x) - (1/n *(sum(Mi$x*ybar_i.$x)))^2)
var_cap_est <- (N^2 * (1/n - 1/N) * sb_2 ) + (N/n * sum(Mi$x^2 * (1/mi - 1/Mi$x) * sw_i_2))
rel_var_y <- var_cap_est/Y_cap_est
hasil <- cbind(Y_cap_est,var_cap_est,rel_var_y)
return(hasil)
}
r102 = unique(dt2$kab)
datakoeh = list()
for(b in r102){
datakoeh[[b]] <- subset(dt2, kab %in% b)
}
try_ye <- list()
for(c in r102){
try_ye[[c]] <- sp_2_ct(data = datakoeh[[c]],f1,f2)
}
this_it = data.frame()
for (q in r102) {
this_it <- rbind(this_it, try_ye[[q]])
}
this_it = cbind(r102, this_it)
}
#fungsi estimasi synthetic dan composite
Y_synth_comp <- function(x,y,vary,Ni){
#synthetic function
N <- Ni
xty <- t(x)%*%y
xtx <- t(x)%*%x
bt_xy <- solve(xtx)%*%xty
bt_xy = t(bt_xy)
Y_syn <- c(rep(0,nrow(x)))
for(i in 1: ncol(x)){
Y_syn = Y_syn + x[,i]*bt_xy[,i]
}
m <- length(Y_syn)
Y_syn <- as.data.frame(Y_syn)
mse <- sum((y - Y_syn)^2)/m-2
var_y_cap_is <- c(rep(mse,m))
Ni_2 <- N^2
ba.1 <- (sum((Y_syn-y)^2/Ni_2))/m
ba.2 <- (sum(rep(var(Y_syn-y),m)/Ni_2))/m
ba_a <- ba.1 - ba.2
mse_syn <- var_y_cap_is + (Ni_2*ba_a)
rel_mse_syn <- mse_syn/Y_syn
#composite function
tetha_i <- mse_syn/(vary+mse_syn)
tetha_i[tetha_i<0]<-0
tetha_i[tetha_i>1]<-1
Y_comp <- (tetha_i*y) + ((1-tetha_i)*Y_syn)
mse_comp <- (1-tetha_i)*mse_syn
rel_mse_comp <- mse_comp/Y_comp
Syn_Comp <- data.frame(cbind(Y_syn, Y_comp,mse_syn, mse_comp, rel_mse_syn, rel_mse_comp))
colnames(Syn_Comp) = c("Y syn","Y comp", "MSE syn", "MSE comp", "Relatif MSE syn","Relatif MSE comp")
return(Syn_Comp)
}
dataset = list()
est_dataset = list()
for( j in 1:iterasi){
set.seed(j)
dataset[[j]] = est_2st_cl(dataY, f1, f2)
est_dataset [[j]] = Y_synth_comp(dataX,dataset[[j]]$Y_cap_est,dataset[[j]]$var_cap_est,Ni)
}
databind = data.frame()
databind_ind = data.frame()
for (b in 1:iterasi) {
databind = rbind(databind, dataset[[b]])
databind_ind = rbind(databind_ind, est_dataset[[b]])
}
datafinal = aggregate(list(databind$Y_cap_est,databind_ind$`Y syn`,databind_ind$`Y comp`,
databind$var_cap_est,databind_ind$`MSE syn`,databind_ind$`MSE comp`
,databind$rel_var_y,databind_ind$`Relatif MSE syn`,databind_ind$`Relatif MSE comp`),
by=list(databind$r102),
FUN=mean)
colnames(datafinal) = c("kab", "Y dir", "Y syn","Y com", "MSE dir","MSE syn","MSE com",
"Relatif MSE dir","Relatif MSE syn","Relatif MSE com")
datafinal$Efisiensi_syn <- datafinal$`MSE syn`/datafinal$`MSE dir`
datafinal$Efisiensi_com <- datafinal$`MSE com`/datafinal$`MSE dir`
return(datafinal)
}Fungsi Estimasi untuk metode Two stage cluster PPSWR-SRSWOR
#Pakai rumus awal
resampling_2st_pps = function(dataY, dataX, muatan, f1, f2, iterasi ){
est_2st_cl <- function(dataY, f1,f2){
sp_2_ct <- function(dataY,f1,f2){
kb_1 <- dataY
kb_1$bs_o <- as.character(paste(kb_1$bs_o))
Mi_kb_1 <- data.frame(table(kb_1$bs_o))
colnames(Mi_kb_1) <- c("bs_o","Muatan_x")
samp_th1 <- function(data,f1){ #pengambilan sampel tahap 1
pps_wr <- function(sizes,n){
N <- length(sizes)
cumsizes <- cumsum(sizes)
totsize <- cumsizes[N]
s <- numeric(n)
for (u in 1:n) {
r <- runif(1,0,totsize)
i <- 1
while (cumsizes[i] < r) {
i <- i+1
}
s[u] <- i
}
return(s)
}
n_length <- f1*nrow(data)
data_sa <- data[pps_wr(data$Muatan_x,n_length),]
return(data_sa)
}
hs_sp_t1 <- samp_th1(Mi_kb_1,f1)
samp_th2 <- function(dt_kb_1,hs_sp_1,f2){
index <- c(1:nrow(hs_sp_1))
hs_sp_1 <- cbind(index,hs_sp_1)
dt_2 <- list();dt_3 <- list() ; samp_rt <- list() ; kode <- list() ; mi <- list()
hs_samp_srs <- list() ; hs_smp_srs <- list(); n_samp <- list()
for(i in hs_sp_1$index){ #Pengambilan sampel ssu scr srs wor
dt_2[[i]] <- subset(hs_sp_1, index %in% i)
dt_3[[i]] <- subset(dt_kb_1, bs_o %in% dt_2[[i]]$bs_o)
samp_rt[[i]] <- sample(1:length(dt_3[[i]]$pop), size = f2*length(dt_3[[i]]$pop), replace = FALSE)
hs_samp_srs[[i]] <- dt_3[[i]][samp_rt[[i]],]
n_samp [[i]] <- nrow(hs_samp_srs[[i]])
kode[[i]] <- rep(i,n_samp[[i]])
mi[[i]] <- rep(n_samp[[i]],n_samp[[i]])
hs_smp_srs[[i]] <- cbind("kode"=kode[[i]],hs_samp_srs[[i]], "mi"=mi[[i]] )
}
gb_samp = data.frame()
for (j in hs_sp_1$index) {
gb_samp <- rbind(gb_samp,hs_smp_srs[[j]])
}
return(gb_samp)
}
res_sp <- samp_th2(kb_1,hs_sp_t1,f2) #kab 1
res_sp$pop <- as.numeric(paste(res_sp$pop))
res_sp$mi <- as.numeric(paste(res_sp$mi))
y_i. <- aggregate(res_sp$pop, list(res_sp$kode), mean)
M_nol <- sum(Mi_kb_1$Muatan_x)
Y_cap_est <- 1/ (length(y_i.$x)) * (sum(y_i.$x * M_nol))
pk_1 <- M_nol*y_i.$x
pk_2 <- pk_1 - Y_cap_est
pk_2_x <- pk_2^2
var_cap_est <- (sum(pk_2_x))/(((length(y_i.$x)))*((length(y_i.$x))-1))
rel_var_y <- var_cap_est/Y_cap_est
hasil <- cbind(Y_cap_est,var_cap_est,rel_var_y)
return(hasil)
}
r102 = unique(dt2$kab)
datakoeh = list()
for(b in r102){
datakoeh[[b]] <- subset(dt2, kab %in% b)
}
try_ye <- list()
for(c in r102){
try_ye[[c]] <- sp_2_ct(data = datakoeh[[c]],f1,f2)
}
this_it = data.frame()
for (q in r102) {
this_it <- rbind(this_it, try_ye[[q]])
}
this_it = cbind(r102, this_it)
}
#fungsi estimasi synthetic dan composite
Y_synth_comp <- function(x,y,vary){
#synthetic function
xty <- t(x)%*%y
xtx <- t(x)%*%x
bt_xy <- solve(xtx)%*%xty
bt_xy = t(bt_xy)
Y_syn <- c(rep(0,nrow(x)))
for(i in 1: ncol(x)){
Y_syn = Y_syn + x[,i]*bt_xy[,i]
}
Y_syn <- as.data.frame(Y_syn)
mse_syn <- ((Y_syn-y)^2)-vary
rel_mse_syn <- mse_syn/Y_syn
#composite function
tetha_i <- mse_syn/(vary+mse_syn)
tetha_i[tetha_i<0]<-0
tetha_i[tetha_i>1]<-1
Y_comp <- (tetha_i*y) + ((1-tetha_i)*Y_syn)
mse_comp <- (1-tetha_i)*mse_syn
rel_mse_comp <- mse_comp/Y_comp
Syn_Comp <- data.frame(cbind(Y_syn, Y_comp,mse_syn, mse_comp, rel_mse_syn, rel_mse_comp))
colnames(Syn_Comp) = c("Y syn","Y comp", "MSE syn", "MSE comp", "Relatif MSE syn","Relatif MSE comp")
return(Syn_Comp)
}
dataset = list()
est_dataset = list()
for( j in 1:iterasi){
set.seed(j)
dataset[[j]] = est_2st_cl(dataY, f1, f2)
est_dataset [[j]] = Y_synth_comp(dataX,dataset[[j]]$Y_cap_est,dataset[[j]]$var_cap_est)
}
databind = data.frame()
databind_ind = data.frame()
for (b in 1:iterasi) {
databind = rbind(databind, dataset[[b]])
databind_ind = rbind(databind_ind, est_dataset[[b]])
}
datafinal = aggregate(list(databind$Y_cap_est,databind_ind$`Y syn`,databind_ind$`Y comp`,
databind$var_cap_est,databind_ind$`MSE syn`,databind_ind$`MSE comp`
,databind$rel_var_y,databind_ind$`Relatif MSE syn`,databind_ind$`Relatif MSE comp`),
by=list(databind$r102),
FUN=mean)
colnames(datafinal) = c("kab", "Y dir", "Y syn","Y com", "MSE dir","MSE syn","MSE com",
"Relatif MSE dir","Relatif MSE syn","Relatif MSE com")
datafinal$Efisiensi_syn <- datafinal$`MSE syn`/datafinal$`MSE dir`
datafinal$Efisiensi_com <- datafinal$`MSE com`/datafinal$`MSE dir`
return(datafinal)
}b. MEnggunakan rumus MSE synthetic yang perbaikanresampling_2st_pps_2 = function(dataY, dataX, Ni, f1, f2, iterasi ){
est_2st_cl <- function(dataY, f1,f2){
sp_2_ct <- function(dataY,f1,f2){
kb_1 <- dataY
kb_1$bs_o <- as.character(paste(kb_1$bs_o))
Mi_kb_1 <- data.frame(table(kb_1$bs_o))
colnames(Mi_kb_1) <- c("bs_o","Muatan_x")
samp_th1 <- function(data,f1){ #pengambilan sampel tahap 1
pps_wr <- function(sizes,n){
N <- length(sizes)
cumsizes <- cumsum(sizes)
totsize <- cumsizes[N]
s <- numeric(n)
for (u in 1:n) {
r <- runif(1,0,totsize)
i <- 1
while (cumsizes[i] < r) {
i <- i+1
}
s[u] <- i
}
return(s)
}
n_length <- f1*nrow(data)
data_sa <- data[pps_wr(data$Muatan_x,n_length),]
return(data_sa)
}
hs_sp_t1 <- samp_th1(Mi_kb_1,f1)
samp_th2 <- function(dt_kb_1,hs_sp_1,f2){
index <- c(1:nrow(hs_sp_1))
hs_sp_1 <- cbind(index,hs_sp_1)
dt_2 <- list();dt_3 <- list() ; samp_rt <- list() ; kode <- list() ; mi <- list()
hs_samp_srs <- list() ; hs_smp_srs <- list(); n_samp <- list()
for(i in hs_sp_1$index){ #Pengambilan sampel ssu scr srs wor
dt_2[[i]] <- subset(hs_sp_1, index %in% i)
dt_3[[i]] <- subset(dt_kb_1, bs_o %in% dt_2[[i]]$bs_o)
samp_rt[[i]] <- sample(1:length(dt_3[[i]]$pop), size = f2*length(dt_3[[i]]$pop), replace = FALSE)
hs_samp_srs[[i]] <- dt_3[[i]][samp_rt[[i]],]
n_samp [[i]] <- nrow(hs_samp_srs[[i]])
kode[[i]] <- rep(i,n_samp[[i]])
mi[[i]] <- rep(n_samp[[i]],n_samp[[i]])
hs_smp_srs[[i]] <- cbind("kode"=kode[[i]],hs_samp_srs[[i]], "mi"=mi[[i]] )
}
gb_samp = data.frame()
for (j in hs_sp_1$index) {
gb_samp <- rbind(gb_samp,hs_smp_srs[[j]])
}
return(gb_samp)
}
res_sp <- samp_th2(kb_1,hs_sp_t1,f2) #kab 1
res_sp$pop <- as.numeric(paste(res_sp$pop))
res_sp$mi <- as.numeric(paste(res_sp$mi))
y_i. <- aggregate(res_sp$pop, list(res_sp$kode), mean)
M_nol <- sum(Mi_kb_1$Muatan_x)
Y_cap_est <- 1/ (length(y_i.$x)) * (sum(y_i.$x * M_nol))
pk_1 <- M_nol*y_i.$x
pk_2 <- pk_1 - Y_cap_est
pk_2_x <- pk_2^2
var_cap_est <- (sum(pk_2_x))/(((length(y_i.$x)))*((length(y_i.$x))-1))
rel_var_y <- var_cap_est/Y_cap_est
hasil <- cbind(Y_cap_est,var_cap_est,rel_var_y)
return(hasil)
}
r102 = unique(dt2$kab)
datakoeh = list()
for(b in r102){
datakoeh[[b]] <- subset(dt2, kab %in% b)
}
try_ye <- list()
for(c in r102){
try_ye[[c]] <- sp_2_ct(data = datakoeh[[c]],f1,f2)
}
this_it = data.frame()
for (q in r102) {
this_it <- rbind(this_it, try_ye[[q]])
}
this_it = cbind(r102, this_it)
}
#fungsi estimasi synthetic dan composite
Y_synth_comp <- function(x,y,vary,Ni){
#synthetic function
N <- Ni
xty <- t(x)%*%y
xtx <- t(x)%*%x
bt_xy <- solve(xtx)%*%xty
bt_xy = t(bt_xy)
Y_syn <- c(rep(0,nrow(x)))
for(i in 1: ncol(x)){
Y_syn = Y_syn + x[,i]*bt_xy[,i]
}
m <- length(Y_syn)
Y_syn <- as.data.frame(Y_syn)
mse <- sum((y - Y_syn)^2)/m-2
var_y_cap_is <- c(rep(mse,m))
Ni_2 <- N^2
ba.1 <- (sum((Y_syn-y)^2/Ni_2))/m
ba.2 <- (sum(rep(var(Y_syn-y),m)/Ni_2))/m
ba_a <- ba.1 - ba.2
mse_syn <- var_y_cap_is + (Ni_2*ba_a)
rel_mse_syn <- mse_syn/Y_syn
#composite function
tetha_i <- mse_syn/(vary+mse_syn)
tetha_i[tetha_i<0]<-0
tetha_i[tetha_i>1]<-1
Y_comp <- (tetha_i*y) + ((1-tetha_i)*Y_syn)
mse_comp <- (1-tetha_i)*mse_syn
rel_mse_comp <- mse_comp/Y_comp
Syn_Comp <- data.frame(cbind(Y_syn, Y_comp,mse_syn, mse_comp, rel_mse_syn, rel_mse_comp))
colnames(Syn_Comp) = c("Y syn","Y comp", "MSE syn", "MSE comp", "Relatif MSE syn","Relatif MSE comp")
return(Syn_Comp)
}
dataset = list()
est_dataset = list()
for( j in 1:iterasi){
set.seed(j)
dataset[[j]] = est_2st_cl(dataY, f1, f2)
est_dataset [[j]] = Y_synth_comp(dataX,dataset[[j]]$Y_cap_est,dataset[[j]]$var_cap_est,Ni)
}
databind = data.frame()
databind_ind = data.frame()
for (b in 1:iterasi) {
databind = rbind(databind, dataset[[b]])
databind_ind = rbind(databind_ind, est_dataset[[b]])
}
datafinal = aggregate(list(databind$Y_cap_est,databind_ind$`Y syn`,databind_ind$`Y comp`,
databind$var_cap_est,databind_ind$`MSE syn`,databind_ind$`MSE comp`
,databind$rel_var_y,databind_ind$`Relatif MSE syn`,databind_ind$`Relatif MSE comp`),
by=list(databind$r102),
FUN=mean)
colnames(datafinal) = c("kab", "Y dir", "Y syn","Y com", "MSE dir","MSE syn","MSE com",
"Relatif MSE dir","Relatif MSE syn","Relatif MSE com")
datafinal$Efisiensi_syn <- datafinal$`MSE syn`/datafinal$`MSE dir`
datafinal$Efisiensi_com <- datafinal$`MSE com`/datafinal$`MSE dir`
return(datafinal)
}