Define Libraries

library("stringr")
library("MASS")
library("dplyr")
library("car")
library("knitr")
library("sandwich")
library("lmtest")
library("bucky")

Define Path

dir.wrk <- getwd()
dir.data <- file.path(dir.wrk, "dataset/data")
dir.annot <- file.path(dir.wrk, "dataset/annotation")

Define Files

file.dat <- file.path(dir.data, "mastertbl_household_reconstruction.tsv")

Load Data

dat <- read.delim(file.dat, header = TRUE, stringsAsFactors = FALSE)

# FACTORIZE ---
dat$DISTRICT <- factor(dat$DISTRICT, levels = c("Dhading", "Sindhupalchowk"))
dat$VDC_MUN <- factor(dat$VDC_MUN, levels = c("Gumdi", "Sertung", "Batase", "Pangtang"))
dat$WARD <- as.factor(dat$WARD)
dat$SEX <- factor(dat$SEX, levels = c("Male", "Female"))
dat$ETHNICITY <- factor(dat$ETHNICITY, levels = c("Bahun", "Chhetri", "Dalit-Others", 
    "Gurung-Magar", "Newar", "Tamang"))
dat$EDUCATION_LEVEL <- factor(dat$EDUCATION_LEVEL, levels = c("Illiterate", "Primary", 
    "Secondary", "University"))
dat$OCCUPATION <- factor(dat$OCCUPATION, levels = c("Agriculture", "Business", "Goverment Employee", 
    "Labour", "Teacher"))
dat$INCOME_LEVEL <- factor(dat$INCOME_LEVEL, levels = c("0-2.5K", "2.5-10K", "10-20K", 
    "20-40K", "40-60K", "60K-Above"))
dat$RECON_COMPLETED <- as.factor(dat$RECON_COMPLETED)
dat$SATISFACTION_LEVEL <- as.factor(dat$SATISFACTION_LEVEL)
dat$FREEDOM_CHOICE <- as.factor(dat$FREEDOM_CHOICE)

# str(dat)



Bootstrapping of Odds ratio with confidence intervals

# FIT ORDERED LOGISTIC REGRESSION ---
fit.m <- MASS::polr(SATISFACTION_LEVEL ~ SEX+ETHNICITY+FREEDOM_CHOICE, data=dat, Hess=TRUE, method="logistic")

# BOOTSTRAPPING
or_boot <- suppressWarnings(car::Boot(fit.m, f = function(fit.m) exp(coef(fit.m)), R=1000))
## Loading required namespace: boot
s_or_boot <- summary(or_boot)

knitr::kable(format(as.data.frame(s_or_boot), digits=3), align="lllll")
R original bootBias bootSE bootMed
SEXFemale 1000 0.680 1.44e-02 1.39e-01 0.680
ETHNICITYChhetri 1000 0.652 4.81e-02 3.06e-01 0.643
ETHNICITYDalit-Others 1000 0.557 3.36e-02 2.37e-01 0.554
ETHNICITYGurung-Magar 1000 0.427 1.04e-02 1.65e-01 0.409
ETHNICITYNewar 1000 0.539 3.07e-02 2.45e-01 0.530
ETHNICITYTamang 1000 0.275 1.38e-03 9.77e-02 0.262
FREEDOM_CHOICE2 1000 2.973 1.72e+04 5.45e+05 3.234
FREEDOM_CHOICE3 1000 2.463 4.85e-01 1.80e+00 2.490
FREEDOM_CHOICE4 1000 2.321 3.61e-01 1.48e+00 2.357
FREEDOM_CHOICE5 1000 2.728 3.99e-01 1.43e+00 2.891
FREEDOM_CHOICE6 1000 6.659 1.17e+00 3.64e+00 7.149
FREEDOM_CHOICE7 1000 9.446 1.86e+00 5.14e+00 10.223
FREEDOM_CHOICE8 1000 14.268 3.66e+00 1.05e+01 15.413
FREEDOM_CHOICE9 1000 39.671 1.49e+01 4.14e+01 43.933
FREEDOM_CHOICE10 1000 29.897 1.06e+01 3.05e+01 31.446 
df.or_boot <- format(as.data.frame(confint(or_boot, type="bca")), digits=3, scientific=FALSE)

knitr::kable(format(df.or_boot, digits=3), align="ll")
2.5 % 97.5 %
SEXFemale 0.471 1.030
ETHNICITYChhetri 0.279 1.623
ETHNICITYDalit-Others 0.247 1.203
ETHNICITYGurung-Magar 0.207 0.889
ETHNICITYNewar 0.244 1.282
ETHNICITYTamang 0.144 0.547
FREEDOM_CHOICE2 0.066 52.605
FREEDOM_CHOICE3 0.684 7.084
FREEDOM_CHOICE4 0.758 6.245
FREEDOM_CHOICE5 0.945 5.795
FREEDOM_CHOICE6 2.080 14.070
FREEDOM_CHOICE7 2.932 20.135
FREEDOM_CHOICE8 4.116 36.423
FREEDOM_CHOICE9 9.272 113.550
FREEDOM_CHOICE10 7.866 109.405

Bootstrapping for Coeffcients with 95% CI (Type= bca)

beta_boot <- suppressWarnings(car::Boot(fit.m, R=1000))
s_beta_boot <- summary(beta_boot)

knitr::kable(format(as.data.frame(s_beta_boot), digits=3), align="lllll")
R original bootBias bootSE bootMed
SEXFemale 1000 -0.385 -0.00248 0.209 -0.391
ETHNICITYChhetri 1000 -0.427 -0.02874 0.397 -0.465
ETHNICITYDalit-Others 1000 -0.585 -0.01731 0.384 -0.596
ETHNICITYGurung-Magar 1000 -0.852 -0.03074 0.356 -0.881
ETHNICITYNewar 1000 -0.617 -0.03875 0.416 -0.668
ETHNICITYTamang 1000 -1.290 -0.04514 0.354 -1.325
FREEDOM_CHOICE2 1000 1.090 -0.06498 1.962 1.062
FREEDOM_CHOICE3 1000 0.901 0.02323 0.558 0.911
FREEDOM_CHOICE4 1000 0.842 0.03965 0.514 0.870
FREEDOM_CHOICE5 1000 1.004 0.04148 0.422 1.030
FREEDOM_CHOICE6 1000 1.896 0.07153 0.413 1.965
FREEDOM_CHOICE7 1000 2.246 0.08251 0.420 2.344
FREEDOM_CHOICE8 1000 2.658 0.10316 0.510 2.793
FREEDOM_CHOICE9 1000 3.681 0.15970 0.580 3.849
FREEDOM_CHOICE10 1000 3.398 0.09948 0.659 3.502 
df.beta_boot <- format(as.data.frame(confint(beta_boot, type="bca")), digits=3, scientific=FALSE)

knitr::kable(format(df.beta_boot, digits=3), align="ll")
2.5 % 97.5 %
SEXFemale -0.767 0.0522
ETHNICITYChhetri -1.125 0.3961
ETHNICITYDalit-Others -1.363 0.1610
ETHNICITYGurung-Magar -1.567 -0.2139
ETHNICITYNewar -1.342 0.2882
ETHNICITYTamang -1.976 -0.5798
FREEDOM_CHOICE2 -2.465 4.1696
FREEDOM_CHOICE3 -0.229 1.9457
FREEDOM_CHOICE4 -0.160 1.8869
FREEDOM_CHOICE5 0.139 1.8419
FREEDOM_CHOICE6 0.985 2.6353
FREEDOM_CHOICE7 1.225 2.9486
FREEDOM_CHOICE8 1.482 3.4657
FREEDOM_CHOICE9 2.356 4.6161
FREEDOM_CHOICE10 2.077 4.6824