Regression Analysis: Happiness and its Potential Predictors
Define Libraries
library("stringr")
library("MASS")
library("dplyr")
library("car")
library("knitr")
library("effects")
library("ggplot2")
library("gridExtra")
library("brant")
library("VGAM")
library("sandwich")
library("lmtest")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)Ordered Logistic Regression between Satisfaction Level and its all potential predictors
# FIT ORDERED LOGISTIC REGRESSION ---
fit.mc <- MASS::polr(SATISFACTION_LEVEL ~ SEX+ETHNICITY+FREEDOM_CHOICE+EDUCATION_LEVEL+OCCUPATION+RECON_COMPLETED+DISTRICT+INCOME_LEVEL+AGE, data=dat, Hess=TRUE, method="logistic")
# CALCULATE P-VALUE ---
summary_table <- coef(summary(fit.mc))
pval <- pnorm(abs(summary_table[, "t value"]), lower.tail=FALSE)* 2
summary_table <- cbind(summary_table, "p value"=round(pval,4))
# CALCULATE LOG-ODDS RATIO ---
tbl_oddratio <- exp(cbind(OddRatio=coef(fit.mc), confint(fit.mc, level=0.95)))
# AGGREGATE DATA ---
summary_table_var <- as.data.frame(cbind(summary_table[1:nrow(tbl_oddratio),], tbl_oddratio))
rmarkdown::paged_table(format(summary_table_var, digits=3), options=list(rows.print=10)) Analysis of deviance
# ANOVA
car::Anova(fit.mc)## Analysis of Deviance Table (Type II tests)
##
## Response: SATISFACTION_LEVEL
## LR Chisq Df Pr(>Chisq)
## SEX 5.036 1 0.024823 *
## ETHNICITY 18.175 5 0.002734 **
## FREEDOM_CHOICE 103.834 9 < 2.2e-16 ***
## EDUCATION_LEVEL 1.097 3 0.777837
## OCCUPATION 3.252 4 0.516514
## RECON_COMPLETED 4.139 4 0.387569
## DISTRICT 2.679 1 0.101701
## INCOME_LEVEL 4.282 5 0.509548
## AGE 0.419 1 0.517282
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1The above analysis reveals that the effects of variables such as, location, education, income, physical level of reconstruction, and occupation on the self-reported level of happiness are not statistically significant (at 95% confidence level). However, the effects of variables such as, Sex, Ethnicity, and Freedom of Choice on the happiness level were found to statistically significant. Thus, we redo the above analysis considering only Sex, Ethinicity, and Freedom of Choice.
Ordered Logistic Regression between Satisfaction Level and its potential predictors
# FIT ORDERED LOGISTIC REGRESSION ---
fit.m <- MASS::polr(SATISFACTION_LEVEL ~ SEX+ETHNICITY+FREEDOM_CHOICE, data=dat, Hess=TRUE, method="logistic")
# CALCULATE P-VALUE ---
summary_table <- coef(summary(fit.m))
pval <- pnorm(abs(summary_table[, "t value"]), lower.tail=FALSE)* 2
summary_table <- cbind(summary_table, "p value"=round(pval,4))
# CALCULATE LOG-ODDS RATIO ---
tbl_oddratio <- exp(cbind(OddRatio=coef(fit.m), confint(fit.m, level=0.95)))
# AGGREGATE DATA ---
summary_table_var <- as.data.frame(cbind(summary_table[1:nrow(tbl_oddratio),], tbl_oddratio))
rmarkdown::paged_table(format(summary_table_var, digits=3), options=list(rows.print=10)) # SATISFACTION THRESHOLD ---
summary_table_st <- as.data.frame(summary_table[(nrow(tbl_oddratio)+1):nrow(summary_table),])
rmarkdown::paged_table(summary_table_st)Analysis of deviance
# ANOVA
car::Anova(fit.m)## Analysis of Deviance Table (Type II tests)
##
## Response: SATISFACTION_LEVEL
## LR Chisq Df Pr(>Chisq)
## SEX 4.008 1 0.045273 *
## ETHNICITY 18.659 5 0.002224 **
## FREEDOM_CHOICE 105.780 9 < 2.2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1Heteroscedasticity Robust stdandard errors and p-values
coef_t <- lmtest::coeftest(fit.m, vcov = sandwich::vcovHC(fit.m, type = "HC0"))
coef_t##
## t test of coefficients:
##
## Estimate Std. Error t value Pr(>|t|)
## SEXFemale -0.38499 0.19202 -2.0049 0.0457280 *
## ETHNICITYChhetri -0.42708 0.37903 -1.1268 0.2606014
## ETHNICITYDalit-Others -0.58535 0.36637 -1.5977 0.1109925
## ETHNICITYGurung-Magar -0.85210 0.34197 -2.4918 0.0131646 *
## ETHNICITYNewar -0.61719 0.39302 -1.5704 0.1172113
## ETHNICITYTamang -1.29017 0.33796 -3.8176 0.0001588 ***
## FREEDOM_CHOICE2 1.08959 2.14033 0.5091 0.6110126
## FREEDOM_CHOICE3 0.90142 0.53259 1.6925 0.0914215 .
## FREEDOM_CHOICE4 0.84210 0.48947 1.7204 0.0862221 .
## FREEDOM_CHOICE5 1.00371 0.40531 2.4764 0.0137342 *
## FREEDOM_CHOICE6 1.89603 0.40113 4.7267 3.289e-06 ***
## FREEDOM_CHOICE7 2.24556 0.40357 5.5643 5.174e-08 ***
## FREEDOM_CHOICE8 2.65805 0.47606 5.5835 4.675e-08 ***
## FREEDOM_CHOICE9 3.68063 0.55919 6.5820 1.658e-10 ***
## FREEDOM_CHOICE10 3.39775 0.61403 5.5336 6.082e-08 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1Odds Ratios with robust confidence intervals
robust_sr <- sqrt(diag( sandwich::vcovHC(fit.m, type = "HC0") ))
y_or <- exp(coef(fit.m))
robust_or <- suppressWarnings(
cbind(coef_t, y_or,
LL = exp(coef(fit.m) - 1.96 * robust_sr),
UL = exp(coef(fit.m) + 1.96 * robust_sr)))
knitr::kable(format(as.data.frame(robust_or), digits=3), align="lllllll")| Estimate | Std. Error | t value | Pr(>|t|) | y_or | LL | UL | |
|---|---|---|---|---|---|---|---|
| SEXFemale | -0.385 | 0.192 | -2.005 | 4.57e-02 | 0.680 | 0.4670 | 0.991 |
| ETHNICITYChhetri | -0.427 | 0.379 | -1.127 | 2.61e-01 | 0.652 | 0.3104 | 1.371 |
| ETHNICITYDalit-Others | -0.585 | 0.366 | -1.598 | 1.11e-01 | 0.557 | 0.2716 | 1.142 |
| ETHNICITYGurung-Magar | -0.852 | 0.342 | -2.492 | 1.32e-02 | 0.427 | 0.2182 | 0.834 |
| ETHNICITYNewar | -0.617 | 0.393 | -1.570 | 1.17e-01 | 0.539 | 0.2497 | 1.165 |
| ETHNICITYTamang | -1.290 | 0.338 | -3.818 | 1.59e-04 | 0.275 | 0.1419 | 0.534 |
| FREEDOM_CHOICE2 | 1.090 | 2.140 | 0.509 | 6.11e-01 | 2.973 | 0.0448 | 197.285 |
| FREEDOM_CHOICE3 | 0.901 | 0.533 | 1.693 | 9.14e-02 | 2.463 | 0.8672 | 6.996 |
| FREEDOM_CHOICE4 | 0.842 | 0.489 | 1.720 | 8.62e-02 | 2.321 | 0.8893 | 6.059 |
| FREEDOM_CHOICE5 | 1.004 | 0.405 | 2.476 | 1.37e-02 | 2.728 | 1.2328 | 6.038 |
| FREEDOM_CHOICE6 | 1.896 | 0.401 | 4.727 | 3.29e-06 | 6.659 | 3.0338 | 14.618 |
| FREEDOM_CHOICE7 | 2.246 | 0.404 | 5.564 | 5.17e-08 | 9.446 | 4.2827 | 20.833 |
| FREEDOM_CHOICE8 | 2.658 | 0.476 | 5.583 | 4.68e-08 | 14.268 | 5.6124 | 36.275 |
| FREEDOM_CHOICE9 | 3.681 | 0.559 | 6.582 | 1.66e-10 | 39.671 | 13.2582 | 118.705 |
| FREEDOM_CHOICE10 | 3.398 | 0.614 | 5.534 | 6.08e-08 | 29.897 | 8.9734 | 99.608 |
Happiness/Satistfaction level for each level of freedom of choice
# EFFECT OF THE FITTED MODEL ---
ef1 <- effects::Effect(focal.predictors="FREEDOM_CHOICE", mod=fit.m,
xlevels=list(FREEDOM_CHOICE = 1:10, 2), latent = TRUE)
#plot(ef1, rug=FALSE)#### FUNCTION: plotEffect()
plotEffect <- function(ef.obj){
# PREPARE DATA ---
df <- cbind(ef.obj$x, fit=ef.obj$fit, upper=ef.obj$upper, lower=ef.obj$lower)
dt <- data.frame(Level=names(ef.obj$thresholds), thresholds=ef.obj$thresholds)
# FACTORIZE ---
df$FREEDOM_CHOICE <- as.numeric(as.character(df$FREEDOM_CHOICE))
# PLOT ---
p <- ggplot(data = df, aes(x=FREEDOM_CHOICE, y=fit)) +
geom_hline(yintercept = dt$thresholds[5:9], color="#D9D9D9", size=1, linetype=2) +
#geom_ribbon(aes(ymin=lower, ymax=upper), fill="#D9D9D9", alpha=0.5) +
geom_line(aes(y = fit), color="#E31A1C", size=1.5) +
geom_errorbar(data=df, aes(x=FREEDOM_CHOICE, ymin=upper, ymax=lower), width=0.2, size=1, color="#969696") +
geom_point(color="#E31A1C", size=3) +
scale_x_continuous(breaks=c(1:10), labels=c(1:10)) +
scale_y_continuous(breaks=dt$thresholds, labels=dt$Level, name="Satisfaction Level",
sec.axis=sec_axis(~.*1,name="Effect", breaks=seq(-2,4,1), labels=seq(-2,4,1))) +
theme(
axis.text.x = element_text(size = 10, color="#000000"),
axis.text.y = element_text(size = 10, color="#000000"),
axis.title = element_text(size = 10, color="#000000"),
plot.title = element_text(size = 10, color="#000000", hjust=0.5),
panel.grid.major = element_blank(),
panel.grid.minor = element_blank(),
axis.ticks = element_line(size=0.4, color="#000000"),
strip.text = element_text(size=10, color="#000000"),
strip.background = element_rect(fill="#FFFFFF", color="#FFFFFF"),
panel.background = element_rect(fill="#FFFFFF", color="#000000"),
legend.text = element_text(size = 10, color="#000000"),
legend.title = element_blank(),
legend.key.size = unit(0.5, "cm"),
legend.position = "none") +
ylab("Effect") +
xlab("Freedom of Choice") +
ggtitle("")
return(p)
}plotEffect(ef.obj=ef1)
Happiness/Satistfaction level for freedom of choice by sex and ethnicity
# EFFECT OF THE FITTED MODEL ---
ef2 <- effects::Effect(focal.predictors=c("SEX","FREEDOM_CHOICE","ETHNICITY"), mod=fit.m,
xlevels=list(FREEDOM_CHOICE = 1:10), latent = TRUE)
#plot(ef2, rug = FALSE, ylim = c(-2.0,4.5))# FUNCTION: plotEffect_MultiVar()
plotEffect_MultiVar <- function(ef.obj){
# PREPARE DATA ---
df <- cbind(ef.obj$x, fit=ef.obj$fit, upper=ef.obj$upper, lower=ef.obj$lower)
dt <- data.frame(Level=names(ef.obj$thresholds), thresholds=ef.obj$thresholds)
# FACTORIZE ---
df$FREEDOM_CHOICE <- as.numeric(as.character(df$FREEDOM_CHOICE))
df$SEX <- factor(df$SEX, levels=c("Male","Female"))
df$ETHNICITY <- factor(df$ETHNICITY, levels=c("Bahun","Chhetri","Dalit-Others","Gurung-Magar","Newar","Tamang"))
# PLOT ---
p <- ggplot(data = df, aes(x=FREEDOM_CHOICE, y=fit)) +
geom_hline(yintercept = dt$thresholds[4:9], color="#D9D9D9", size=1, linetype=2) +
#geom_ribbon(aes(ymin=lower, ymax=upper), fill="#D9D9D9", alpha=0.5) +
geom_line(aes(y = fit), color="#E31A1C", size=2) +
geom_errorbar(data=df, aes(x=FREEDOM_CHOICE, ymin=upper, ymax=lower), width=0.2, size=1, color="#969696") +
geom_point(color="#E31A1C", size=2) +
scale_x_continuous(breaks=c(1:10), labels=c(1:10)) +
scale_y_continuous(breaks=dt$thresholds, labels=dt$Level, name="Satisfaction Level",
sec.axis=sec_axis(~.*1,name="Effect", breaks=seq(-2,4,1), labels=seq(-2,4,1))) +
#facet_grid(ETHNICITY~SEX, scales="fixed", space="fixed") +
facet_grid(SEX~ETHNICITY, scales="fixed", space="fixed") +
theme(
axis.text.x = element_text(size = 10, color="#000000"),
axis.text.y = element_text(size = 10, color="#000000"),
axis.title = element_text(size = 10, color="#000000"),
plot.title = element_text(size = 10, color="#000000", hjust=0.5),
panel.grid.major = element_blank(),
panel.grid.minor = element_blank(),
axis.ticks = element_line(size=0.4, color="#000000"),
strip.text = element_text(size=10, color="#000000"),
strip.background = element_rect(fill="#FFFFFF", color="#FFFFFF"),
panel.background = element_rect(fill="#FFFFFF", color="#000000"),
legend.text = element_text(size = 10, color="#000000"),
legend.title = element_blank(),
legend.key.size = unit(0.5, "cm"),
legend.position = "none") +
ylab("Effect") +
xlab("Freedom of Choice") +
ggtitle("")
return(p)
}plotEffect_MultiVar(ef.obj=ef2)
ef3 <- effects::Effect(focal.predictors="FREEDOM_CHOICE", mod=fit.m,
xlevels=list(FREEDOM_CHOICE = 1:10, 2))
#plot(ef3, lines=list(multiline=TRUE), rug=FALSE)Propbability of each level of happiness/satisfaction at different levels of Freedom
#### FUNCTION: plotEffect()
plotProbEffect <- function(ef.obj){
# PREPARE DATA ---
df <- cbind(ef.obj$x, ef.obj$prob)
# RESHAPE DATA ---
dm <- reshape2::melt(df, id.vars="FREEDOM_CHOICE", value.name="PROBABILITY")
colnames(dm) <- c("FREEDOM_CHOICE","SATISFACTION_LEVEL","PROBABILITY")
dm$SATISFACTION_LEVEL <- str_replace_all(dm$SATISFACTION_LEVEL, "prob.X", "")
# FACTORIZE ---
dm$FREEDOM_CHOICE <- as.numeric(as.character(dm$FREEDOM_CHOICE))
dm$SATISFACTION_LEVEL <- as.factor(as.numeric(dm$SATISFACTION_LEVEL))
# COLOR PALETTE
cpalette <- rev(c("#e31a1c","#a6cee3","#1f78b4","#b2df8a","#33a02c",
"#fb9a99","#fdbf6f","#ff7f00","#cab2d6","#6a3d9a"))
# PLOT ---
p <- ggplot(data = dm, aes(x=FREEDOM_CHOICE, y=PROBABILITY)) +
geom_line(aes(y=PROBABILITY, color=SATISFACTION_LEVEL), size=1.5) +
geom_point(aes(y=PROBABILITY, color=SATISFACTION_LEVEL), size=3) +
scale_x_continuous(breaks=c(1:10), labels=c(1:10)) +
scale_y_continuous(breaks=seq(0,0.6,0.1)) +
scale_color_manual(values=cpalette, name="Satisfaction Level") +
theme(
axis.text.x = element_text(size = 10, color="#000000"),
axis.text.y = element_text(size = 10, color="#000000"),
axis.title = element_text(size = 10, color="#000000"),
plot.title = element_text(size = 10, color="#000000", hjust=0.5),
panel.grid.major = element_blank(),
panel.grid.minor = element_blank(),
axis.ticks = element_line(size=0.4, color="#000000"),
strip.text = element_text(size=10, color="#000000"),
strip.background = element_rect(fill="#FFFFFF", color="#FFFFFF"),
panel.background = element_rect(fill="#FFFFFF", color="#000000"),
legend.text = element_text(size = 10, color="#000000"),
legend.title = element_text(size = 10, color="#000000"),
legend.key.size = unit(0.5, "cm"),
legend.background = element_blank(),
legend.box.background = element_blank(),
legend.position = "bottom") +
ylab("Probability of Satisfaction") +
xlab("Freedom of Choice") +
ggtitle("")
return(p)
}plotProbEffect(ef.obj=ef3)
Robustness Check for proportional Odds assumption
suppressWarnings(brant(fit.m))## ----------------------------------------------------
## Test for X2 df probability
## ----------------------------------------------------
## Omnibus 5928.51 120 0
## SEXFemale 27.11 8 0
## ETHNICITYChhetri 2.2 8 0.97
## ETHNICITYDalit-Others 2.17 8 0.98
## ETHNICITYGurung-Magar 2.86 8 0.94
## ETHNICITYNewar 5.42 8 0.71
## ETHNICITYTamang 1.9 8 0.98
## FREEDOM_CHOICE2 19.35 8 0.01
## FREEDOM_CHOICE3 9.99 8 0.27
## FREEDOM_CHOICE4 10.65 8 0.22
## FREEDOM_CHOICE5 10.2 8 0.25
## FREEDOM_CHOICE6 27.36 8 0
## FREEDOM_CHOICE7 25.44 8 0
## FREEDOM_CHOICE8 13.19 8 0.11
## FREEDOM_CHOICE9 4.92 8 0.77
## FREEDOM_CHOICE10 12.01 8 0.15
## ----------------------------------------------------
##
## H0: Parallel Regression Assumption holdsRobustness Check with an Alternative Model: VGAM
m_glm <- suppressWarnings(VGAM::vglm(SATISFACTION_LEVEL~FREEDOM_CHOICE+ETHNICITY+SEX, data = dat,
family = VGAM::cumulative(link = "logitlink", parallel = TRUE, reverse = FALSE)))Anova(m_glm)## Analysis of Deviance Table (Type II tests)
##
## Response: SATISFACTION_LEVEL
## Df Chisq Pr(>Chisq)
## FREEDOM_CHOICE 9 104.8817 < 2.2e-16 ***
## ETHNICITY 5 18.5139 0.002367 **
## SEX 1 3.9999 0.045502 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1summary(m_glm)##
## Call:
## VGAM::vglm(formula = SATISFACTION_LEVEL ~ FREEDOM_CHOICE + ETHNICITY +
## SEX, family = VGAM::cumulative(link = "logitlink", parallel = TRUE,
## reverse = FALSE), data = dat)
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept):1 -4.354025 0.638561 -6.818 9.20e-12 ***
## (Intercept):2 -3.637775 0.533109 -6.824 8.87e-12 ***
## (Intercept):3 -3.125420 0.487341 -6.413 1.42e-10 ***
## (Intercept):4 -2.276158 0.446475 -5.098 3.43e-07 ***
## (Intercept):5 -0.008954 0.431941 -0.021 0.983461
## (Intercept):6 0.742839 0.434555 1.709 0.087372 .
## (Intercept):7 1.517922 0.438535 3.461 0.000537 ***
## (Intercept):8 2.396974 0.445326 5.383 7.35e-08 ***
## (Intercept):9 2.914447 0.451573 6.454 1.09e-10 ***
## FREEDOM_CHOICE2 -1.090869 0.711163 -1.534 0.125049
## FREEDOM_CHOICE3 -0.901337 0.514888 -1.751 0.080024 .
## FREEDOM_CHOICE4 -0.842026 0.426266 -1.975 0.048228 *
## FREEDOM_CHOICE5 -1.003626 0.384186 -2.612 0.008992 **
## FREEDOM_CHOICE6 -1.895944 0.405236 -4.679 2.89e-06 ***
## FREEDOM_CHOICE7 -2.245502 0.407632 -5.509 3.62e-08 ***
## FREEDOM_CHOICE8 -2.657949 0.516012 -5.151 2.59e-07 ***
## FREEDOM_CHOICE9 -3.680532 0.533739 -6.896 5.36e-12 ***
## FREEDOM_CHOICE10 -3.397616 0.432264 -7.860 3.84e-15 ***
## ETHNICITYChhetri 0.427046 0.386501 1.105 0.269202
## ETHNICITYDalit-Others 0.585372 0.393843 1.486 0.137198
## ETHNICITYGurung-Magar 0.852089 0.352033 2.420 0.015500 *
## ETHNICITYNewar 0.617323 0.418800 1.474 0.140474
## ETHNICITYTamang 1.290198 0.339337 3.802 0.000143 ***
## SEXFemale 0.385026 0.192515 2.000 0.045502 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Number of linear predictors: 9
##
## Names of linear predictors: logitlink(P[Y<=1]), logitlink(P[Y<=2]),
## logitlink(P[Y<=3]), logitlink(P[Y<=4]), logitlink(P[Y<=5]), logitlink(P[Y<=6]),
## logitlink(P[Y<=7]), logitlink(P[Y<=8]), logitlink(P[Y<=9])
##
## Residual deviance: 1339.561 on 3405 degrees of freedom
##
## Log-likelihood: -669.7804 on 3405 degrees of freedom
##
## Number of Fisher scoring iterations: 9
##
## Warning: Hauck-Donner effect detected in the following estimate(s):
## '(Intercept):1'
##
##
## Exponentiated coefficients:
## FREEDOM_CHOICE2 FREEDOM_CHOICE3 FREEDOM_CHOICE4
## 0.33592444 0.40602629 0.43083678
## FREEDOM_CHOICE5 FREEDOM_CHOICE6 FREEDOM_CHOICE7
## 0.36654797 0.15017654 0.10587437
## FREEDOM_CHOICE8 FREEDOM_CHOICE9 FREEDOM_CHOICE10
## 0.07009183 0.02520957 0.03345291
## ETHNICITYChhetri ETHNICITYDalit-Others ETHNICITYGurung-Magar
## 1.53272281 1.79565969 2.34454039
## ETHNICITYNewar ETHNICITYTamang SEXFemale
## 1.85395840 3.63350533 1.46965287is.parallel(m_glm)## (Intercept) FREEDOM_CHOICE ETHNICITY SEX
## FALSE TRUE TRUE TRUERobustness Check: OLS Model
dat$SATISFACTION_LEVEL <- as.numeric(dat$SATISFACTION_LEVEL)
fit_lin <- lm(SATISFACTION_LEVEL~SEX+ETHNICITY+FREEDOM_CHOICE, data=dat)
lmtest::coeftest(fit_lin, vcov = sandwich::vcovHC(fit_lin, type = "HC0"))##
## t test of coefficients:
##
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 6.22825 0.43780 14.2264 < 2.2e-16 ***
## SEXFemale -0.24011 0.17847 -1.3454 0.1793355
## ETHNICITYChhetri -0.44600 0.36866 -1.2098 0.2271474
## ETHNICITYDalit-Others -0.73752 0.36165 -2.0393 0.0421384 *
## ETHNICITYGurung-Magar -0.87009 0.34342 -2.5336 0.0117070 *
## ETHNICITYNewar -0.60366 0.38294 -1.5764 0.1158059
## ETHNICITYTamang -1.35496 0.32976 -4.1090 4.911e-05 ***
## FREEDOM_CHOICE2 0.79425 1.13869 0.6975 0.4859261
## FREEDOM_CHOICE3 0.50397 0.49797 1.0120 0.3121921
## FREEDOM_CHOICE4 0.72167 0.43444 1.6611 0.0975450 .
## FREEDOM_CHOICE5 0.86855 0.41681 2.0838 0.0378716 *
## FREEDOM_CHOICE6 1.52355 0.40522 3.7599 0.0001979 ***
## FREEDOM_CHOICE7 1.97191 0.41256 4.7797 2.549e-06 ***
## FREEDOM_CHOICE8 2.34964 0.50468 4.6557 4.527e-06 ***
## FREEDOM_CHOICE9 3.35176 0.48352 6.9321 1.890e-11 ***
## FREEDOM_CHOICE10 2.79334 0.48732 5.7321 2.080e-08 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1