library(foreign)
library(tidyverse)
library(ggplot2)
library(psych)
library(stringr)
library(lme4)
library(performance)
library(car)
library(ggeffects)
library(ggplot2)
library(rsample)
library(glue)
library(caret)
library(randomForest)
library(pROC)
library(mi)
library(doParallel)
library(missForest)
library(VIM)
library(xgboost)
library(SHAPforxgboost)
library(lme4)
library(performance)
library(psych)
library(missMDA)
library(factoextra)
library(sjPlot)
library(forcats)The analysis is organised into three chapters, presented as tabs below. Click any tab to jump to that chapter. Within each chapter, the left-hand table of contents lets you navigate sub-sections. Code blocks are folded by default — click ▶ Code to expand them.
This section invlove cleaning steps which do not involve information on the whole dataset therefore are empty of data leakage. Here we mainly change variable names, group categories, and changes variables types to fit our purpose.
First we load the Afrobarometer, and change the variable names for a better leasibility. We select the sample of interest. Respondent who didn’t receive the vaccine, and which respond to whether or not they wish in the future to be vaccinate. Respondents already vaccinated were “not applicable” for this question, therefore are filtered out. Respondent who were not vaccinate and refuse to answer (a low frequency), are filtered out.
For predictors, we use variables which were collected in all the countries. Some variable such as trust in parliament were excluded as they were not collected in some countries.
afro <- read.spss('/Users/PaulMUTAMBA/Documents/Old Pc/KUL QASS/SEM/Assignment/afrobarometer_r9.sav',
use.value.labels = TRUE, to.data.frame = TRUE) %>% select("hesitancy" = Q58B, "lang_home" = Q2, "lang_itrvw" = Q102,
"diff_schl" = Q40B, "feel_rspect_schl" = Q40D,
"diff_hsptl" = Q41B, "feel_rspect_hsptl" = Q41D,
"diff_id" = Q42B, "feel_rspect_id" = Q42D,
"diff_pol" = Q43B, "trst_pres" = Q37A,
"trst_police" = Q37G,"trst_courts" = Q37I, "ctry" = COUNTRY, "age" = Q1,
"gndr" = Q100, "attnd_cmty_meeting" = Q10A,
"raise_issue" = Q10B, "protest_march" = Q10C,
"urban_rural" = URBRUR, "lvl_edu" = Q94,
"religion" = Q95, "trst_relative" = Q86B,
"trst_neighbours" = Q86C, "trst_people_other_religion" = Q86E, "trst_other_people_know" = Q86D,
"enough_food" = Q6A, "enough_clean_water" = Q6B, "enough_medecines" = Q6C,
"enough_fuel_cook" = Q6D, "enough_cash" = Q6E
)
afro1 <- afro %>% filter(!hesitancy %in% c("Not applicable", "Refused", "Don't know")) %>% droplevels()
##----All Categories As Underscore------##
afro_clean <- afro1 %>% mutate(across(everything(), ~ str_to_lower(.) %>% str_trim()))We end with 22598 respondents. As one can see below there is no NA. The for each respondent who refused to answer, the afrobarometer as a response variable — “refused answer”.
print(nrow(afro_clean))[1] 22598sapply(afro_clean, function(x) sum(is.na(x))) hesitancy lang_home
0 0
lang_itrvw diff_schl
0 0
feel_rspect_schl diff_hsptl
0 0
feel_rspect_hsptl diff_id
0 0
feel_rspect_id diff_pol
0 0
trst_pres trst_police
0 0
trst_courts ctry
0 0
age gndr
0 0
attnd_cmty_meeting raise_issue
0 0
protest_march urban_rural
0 0
lvl_edu religion
0 0
trst_relative trst_neighbours
0 0
trst_people_other_religion trst_other_people_know
0 0
enough_food enough_clean_water
0 0
enough_medecines enough_fuel_cook
0 0
enough_cash
0 Next we start the diverse data cleaning steps. We do it by categories or by group of categories.
Here, we clean the various categories of categorical variables before converting them into numerical variables. This ensures a smooth transition from the categorical format to the numerical format without causing any errors.
We group peri-urban and urban together, so that we end-up only with two categories, rural and urban. The choice was mainly due to the low frequency of respondent from peri-urban area, which impacted the significance of the variable.
table(afro_clean$urban_rural)
peri-urban rural urban
100 12124 10374 afro_clean$urban_rural <- ifelse(afro_clean$urban_rural == "peri-urban", "urban", afro_clean$urban_rural)Next we create the Lived Poverty Index based on the Afrobarometer application — “To measure lived poverty, Afrobarometer asks respondents: Over the past year, how often, if ever, have you or anyone in your family gone without: Enough food to eat? Enough clean water for home use? Medicines or medical treatment? Enough fuel to cook your food? A cash income?”
The LPI being a numeric values, we had to transform those variables into numeric. “No answer” where replaced into NAs. The LPI is calculate as the mean across the 5 variables.
##-----LPI Creation--------##
afro_clean <- afro_clean %>%
mutate(across(
starts_with("enough"),
~as.numeric(factor(.x, levels = c("never", "just once or twice", "several times",
"many times", "always")))
))
afro_clean$lpi <- rowMeans(afro_clean %>% select(starts_with("enough")), na.rm = TRUE)
summary(afro_clean$lpi) Min. 1st Qu. Median Mean 3rd Qu. Max. NA's
1.000 1.800 2.600 2.529 3.200 5.000 1 For these variables we grouped respondent who had no contact with any of those institutions, and those who didn’t answered or didn’t know onto a same category—“no contact/answer”. At the same time, for variables related to trust, we group respondents who refused to answer and those who didn’t know into “no answer”.
afro_merged <- afro_clean %>%
mutate(across(c(matches("^(diff_|feel_)"), starts_with("trst_")), as.factor)) %>%
mutate(across(matches("^(diff_|feel_)"), ~ fct_collapse(.,
"no answer" = c("refused", "don't know")
))) %>%
mutate(across(starts_with("trst_"), ~ fct_collapse(.,
"no answer" = c("refused",
"don’t know/haven’t heard enough to say",
"don't know")
)))
##------Change no contact/answer Categories For difficulty and feel respected for services------##
afro_merged <- afro_merged %>%
mutate(across(
starts_with("feel") | starts_with("diff"),
~ as.factor(.) %>%
fct_collapse("no contact/answer" = c("no answer", "not applicable", "no contact", "don’t know"))
))Here we create four categories for variables related to civic engagement. Active:Combine all “Yes” frequencies, Willing: Keep the “Would do if had chance” group, Resistant: Keep the “Never” group, No Answer: Combine Refused/Don’t know.
afro_merged <- afro_merged %>% # Or afro_edu or afro_relig (whichever is your latest)
mutate(across(c(attnd_cmty_meeting, raise_issue, protest_march), ~ fct_collapse(.,
"Active" = c("yes, once or twice",
"yes, several times",
"yes, often"),
"Willing" = c("no, would do it if i had the chance"),
"Resistant" = c("no, would never do this"),
"No Answer" = c("refused", "don't know")
)))afro_merged <- afro_merged %>% mutate(across(
starts_with("feel"),
~ factor(., levels = c("not at all", "a little bit", "somewhat", "a lot", "no contact/answer") )
)) %>%
mutate(across(
starts_with("diff"),
~ factor(., levels = c("very easy", "easy", "difficult", "very difficult", "no contact/answer"))
)) %>% mutate(across(
c(attnd_cmty_meeting,raise_issue, protest_march),
~ factor(., levels = c("Active", "Willing", "Resistant", "No Answer"))
))We change into numeric, and respondent who refused to answer are marked as NAs.
afro_merged$age <- ifelse(afro_merged$age %in% c("don't know", "refused"), NA, afro_merged$age)
afro_merged$age <- as.numeric(as.character(afro_merged$age))Here we grouped the languages into two categories: colonial and indigenous/others. Colonial language includes—English, French, Portugese, and Afrikaans. Indigeneous/others languages includes all the myriad african languages mentioned in the Afrobarometer plus Arabic, Swahili, and the different Creoles.
afro_merged$lang_home1 <- ifelse(afro_merged$lang_home %in% c("french", "english", "portuguese", "afrikaans", "afrikaans_duplicated_700"), "colonial", "indigeneous/others")
afro_merged$lang_itrvw1 <- ifelse(afro_merged$lang_itrvw %in% c("french", "english", "portuguese", "afrikaans", "afrikaans_duplicated_700"), "colonial", "indigeneous/others")We grouped the different level of education into five categories: “None or Informal”, “Primary”, “Secondary”, “Higher Education”, and “No answer”.
afro_merged <- afro_merged %>%
mutate(lvl_edu_clean = fct_collapse(lvl_edu,
"None or Informal" = c("no formal schooling",
"informal schooling only (including koranic schooling)"),
"Primary" = c("some primary schooling",
"primary school completed"),
"Secondary" = c("intermediate school or some secondary school / high school",
"secondary school / high school completed"),
"Higher Education" = c("post-secondary qualifications, other than university e.g. a diploma or degree from a polytechnic or college",
"some university",
"university completed",
"post-graduate"),
"No Answer" = c("refused", "don’t know")
)) %>% mutate(lvl_edu_clean = factor(lvl_edu_clean, levels = c("None or Informal",
"Primary",
"Secondary",
"Higher Education",
"No Answer")))We grouped the different religion into five categories: Muslim, Catholic, Protestan/Christian, None/Traditional, and Other/No answer.
afro_merged <- afro_merged %>%
mutate(religion_clean = fct_collapse(religion,
"Muslim" = c("muslim only (i.e., respondents says only “muslim”, without identifying a specific sub-group)",
"sunni only (i.e., respondents says only “sunni” or “sunni muslim”, without identifying a specific sub-group)",
"mouridiya brotherhood",
"tijaniya brotherhood",
"qadiriya brotherhood",
"shia",
"ismaeli",
"ançardine"),
"Catholic" = c("roman catholic"),
"Protestant/Christian" = c("christian only (i.e., respondents says only “christian”, without identifying a specific sub-group)",
"pentecostal (e.g., “born again” and/or “saved”)",
"evangelical",
"anglican",
"lutheran",
"seventh day adventist",
"baptist",
"methodist",
"presbyterian",
"calvinist",
"independent (e.g., “african independent church”)",
"zionist christian church",
"church of christ",
"new apostolic church",
"united church of zambia or ucz",
"orthodox",
"coptic",
"jehovah’s witness",
"mormon",
"christian mission in many lands (cmml)",
"dutch reformed",
"eglise du christianisme céleste",
"faith of unity",
"fifohazana",
"mennonite",
"morovian",
"salvation army",
"salvation army_duplicated_823",
"quaker / friends"),
"None/Traditional" = c("none",
"traditional / ethnic religion",
"atheist (do not believe in a god)",
"agnostic (do not know if there is a god)"),
"Other/No Answer" = c("other",
"refused",
"don't know",
"not asked in the country",
"hindu",
"bahai",
"jewish")
))We collapse the willingness to get vaccinate into two categories—hesitant and willing.
afro_merged$hesitancy <- ifelse(afro_merged$hesitancy %in% c("very unlikely", "somewhat unlikely"), "hesitant", "willing")
afro_merged$hesitancy <- factor(afro_merged$hesitancy, levels = c("hesitant", "willing"))
afro_merged$hesitancy <- relevel(as.factor(afro_merged$hesitancy), ref = "willing")
table(afro_merged$hesitancy)
willing hesitant
9820 12778 We turn all trust variable into numeric with “no answer” as NAs.
afro_merged<- afro_merged %>% mutate(across(starts_with("trst"),
~as.numeric(factor(.x,levels = c("not at all","just a little", "somewhat", "a lot")))))Finally, depending on the degree of integration of indigenous languages into each country’s school system, we classified countries into four categories: limited integration, moderate integration, extensive integration and no integration.
afro_gini2 <- afro_merged %>%
mutate(
indigenous_lang_integration = case_when(
ctry %in% c("botswana", "seychelles", "zimbabwe", "uganda", "south africa",
"senegal", "nigeria", "niger", "namibia", "mauritius",
"kenya", "ghana") ~ "limited integration",
ctry %in% c("zambia", "malawi", "lesotho", "gambia", "ethiopia",
"eswatini", "mozambique", "madagascar", "burkina faso") ~ "moderate integration",
ctry == "tanzania" ~ "extensive integration",
TRUE ~ "no integration"
),
indigenous_lang_integration = factor(indigenous_lang_integration,
levels = c("no integration", "limited integration",
"moderate integration", "extensive integration"))
)
df_ml <- afro_gini2 %>% select("lang_home"= lang_home1,"lang_interview" = lang_itrvw1, diff_schl, feel_rspect_schl, diff_hsptl, feel_rspect_hsptl, diff_id,
feel_rspect_id, diff_pol, trst_pres, trst_police, trst_courts, trst_relative, trst_neighbours, trst_people_other_religion,
trst_other_people_know,attnd_cmty_meeting, raise_issue, protest_march, gndr, age, urban_rural, "lvl_edu" = lvl_edu_clean,
"religion" = religion_clean, lpi, indigenous_lang_integration, hesitancy, ctry)
names(df_ml)[c(3, 5, 7, 9)] <- c("difficulty_schl_services", "difficulty_hsptl_services",
"difficulty_id_services", "difficulty_police_services")
#saveRDS(df_ml, "/Users/PaulMUTAMBA/Documents/Old Pc/KUL QASS/Master/R_data/df_ml.rds")The df_ml data table is saved to prevent any data leaks that might occur during advanced data processing. Consequently, the df_ml data table will be called separately during machine learning and multi-level advanced data processing.
Here we call the df_ml dataset saved during the previous stage. All the Machine Learning Process are done with a seed set to 42.
df <- readRDS('/Users/PaulMUTAMBA/Documents/Old Pc/KUL QASS/Master/R_data/df_ml.rds') %>% ungroup() %>% select(-gdp, -gini)
set.seed(42)Before any preprocessing we split the DataFrame into train and test sets. All the preprocess steps will be applied separately for the train and test set.
split <- initial_split(df, prop = .8)
train_set <- training(split)
test_set <- testing(split)
#saveRDS(test_set, "/Users/PaulMUTAMBA/Documents/Old Pc/KUL QASS/Master/R_data/test_set.rds")To avoid Data Leakage, we saved the test_set separately.
We have employed advanced data processing techniques here. We have converted categorical variables into numerical variables and combined various scores into a single score. This mainly applies to indices relating to institutions and interpersonal relationships, as well as scores for civic engagement. We too collapse the language policy from four to three, as Tanzania is the only country with extensive integration.
diff_cols <- c("difficulty_schl_services", "difficulty_hsptl_services",
"difficulty_id_services", "difficulty_police_services")
respect_cols <- c("feel_rspect_schl", "feel_rspect_hsptl", "feel_rspect_id")
df <- train_set %>%
mutate(
# 1. Overall Service Utilization (1 = Yes, 0 = No)
# Checks if any of the columns are NOT "no contact/answer"
service_utilization = as.integer(if_any(all_of(diff_cols), ~ . != "no contact/answer" & !is.na(.))),
# 2. Convert text to numbers for Difficulty (creates temporary columns)
across(all_of(diff_cols),
~ case_when(
. == "very easy" ~ 1,
. == "easy" ~ 2,
. == "difficult" ~ 3,
. == "very difficult" ~ 4,
TRUE ~ NA_real_ # Everything else, including "no contact", becomes NA
), .names = "{.col}_num"),
# 3. Convert text to numbers for Respect (creates temporary columns)
across(all_of(respect_cols),
~ case_when(
. == "not at all" ~ 1,
. == "a little bit" ~ 2,
. == "somewhat" ~ 3,
. == "a lot" ~ 4,
TRUE ~ NA_real_
), .names = "{.col}_num")
) %>%
# 4. Calculate the row means, ignoring NAs
rowwise() %>%
mutate(
avg_difficulty = mean(c_across(ends_with("_num") & contains("difficulty")), na.rm = TRUE),
avg_respect = mean(c_across(ends_with("_num") & contains("rspect")), na.rm = TRUE)
) %>%
ungroup() %>%
# Drop the temporary numeric columns to keep the dataframe clean
select(-ends_with("_num"), -starts_with("diff"), -starts_with("feel"))
df <- df %>%
mutate(
avg_diff_model = if_else(is.na(avg_difficulty), 0, avg_difficulty),
avg_resp_model = if_else(is.na(avg_respect), 0, avg_respect)
) %>% select(-avg_difficulty, -avg_respect)
##------Trust Transformation-------##
df_mlvl <- df
##----Creation of country level variables have information on the whole dataset----------##
df_mlvl$service_utilization <- as.factor(ifelse(df_mlvl$service_utilization == 1, "yes", "no"))
#df_mlvl <- df_mlvl %>% mutate(across(where(is.numeric), ~as.numeric(scale(.))))
df_mlvl$service_utilization <- ifelse(df_mlvl$service_utilization == "yes", 1, 0)
col_num <- colnames(df_mlvl %>% select_if(~class(.) == "numeric"))
df_mlvl$Interpersonal_Index <- rowMeans(df_mlvl[, c("trst_relative", "trst_neighbours", "trst_people_other_religion", "trst_other_people_know")], na.rm=TRUE)
df_mlvl$Institutional_Index <- rowMeans(df_mlvl[, c("trst_pres", "trst_police", "trst_courts")], na.rm=TRUE)
df_mlvl <- df_mlvl %>%
mutate(
language_profile = case_when(
lang_home == "colonial" & lang_interview == "colonial" ~ "Consistent Colonial",
lang_home == "indigeneous/others" & lang_interview == "indigeneous/others" ~ "Consistent Indigenous",
(lang_home == "colonial" & lang_interview == "indigeneous/others") |
(lang_home == "indigeneous/others" & lang_interview == "colonial") ~ "Mixed/Code-Switcher",
TRUE ~ NA_character_
)
)
civic_cols <- c("attnd_cmty_meeting", "raise_issue", "protest_march")
df_mlvl <- df_mlvl %>%
mutate(
# 1. Convert text to numbers using case_when (creates temporary columns)
across(all_of(civic_cols),
~ case_when(
. == "Active" ~ 3,
. == "Willing" ~ 2,
. == "Resistant" ~ 1,
TRUE ~ NA_real_ # "No Answer" and anything else becomes NA
), .names = "{.col}_num")
) %>%
# 2. Calculate the row means, ignoring NAs
rowwise() %>%
mutate(
civic_participation_score = mean(c_across(ends_with("_num") & contains(c("attnd", "raise", "protest"))), na.rm = TRUE)
) %>%
ungroup() %>%
# 3. Clean up by dropping the temporary numeric columns
select(-ends_with("_num"))
df_mlvl$lang_policy <- fct_collapse(
df_mlvl$indigenous_lang_integration,
"no integration" = "no integration",
"limited integration" = "limited integration",
"moderate or extensive" = c("moderate integration", "extensive integration")
)
#####_______FIN______###
train_set <- df_mlvl %>% mutate(across(where(~class(.) == "character"), ~as.factor(.))) %>%
select(hesitancy, language_profile, avg_diff_model, avg_resp_model, service_utilization,
Interpersonal_Index, Institutional_Index, civic_participation_score, gndr, age, urban_rural, lvl_edu, religion,
lpi,ctry, lang_policy)train_set is the final set after all the prepossessing steps. The next step consist of Imputing NA’s.
##----imputation------##
unique_countries <- unique(train_set$ctry)
train_set$ctry <- as.factor(train_set$ctry)
country_dfs <- list()
for (country in unique_countries) {
country_df <- train_set[train_set$ctry == country, ]
country_dfs[[country]] <- as.data.frame(country_df)
}
library(mi)
library(doParallel)
library(missForest)
#summary(train_set$avg_diff_model)
fact_var <- colnames(train_set %>% select_if(~class(.) %in% c("character", "factor")))
num_var <- colnames(train_set %>% select_if(~!class(.) %in% c("character", "factor")))
mdf <- list()
for (country in unique_countries) {
mdf[[country]] <- missing_data.frame (as.data.frame(country_dfs[[country]]) )
mdf[[country]] <- change (mdf[[country]] , y = num_var, what = "type", to = "continuous")
mdf[[country]] <- change (mdf[[country]] , y = fact_var, what = "type", to = "ordered-categorical")
}
##-----Randomness--------##
#image(mdf[[2]])
library(VIM)
#country_dfs[[2]] %>% mosaicMiss(highlight = "trst_courts" , plotvars = c("religion", "gndr"))
##Conclusion not MCAR, but MAR
##-----MissForest Imputation----##
registerDoParallel(cores=4)
set.seed(42)
imputed_missForests <- list()
for (country in unique_countries){
imputed_missForests[[country]] <- missForest(country_dfs[[country]], verbose = TRUE, parallelize = "forests")
}
imputed_missForests[[1]]$OOBerror[1]
missForest_imputed <- data.frame(
original = country_dfs[[1]]$trst_courts,
imputed_missForest = imputed_missForests[[1]]$ximp$trst_courts
)
train_set_imp <- data.frame()
for (country in unique_countries){
train_set_imp <- rbind(train_set_imp, as.data.frame(imputed_missForests[[country]]$ximp))
}
train_set_imp <- train_set_imp %>% mutate(across(where(is.numeric), round))
#saveRDS(train_set_imp, '/Users/PaulMUTAMBA/Documents/Old Pc/KUL QASS/Master/R_data/train_set_imp2.rds')The train_set_imp is the final set after imputation.
train_set_imp <- readRDS('/Users/PaulMUTAMBA/Documents/Old Pc/KUL QASS/Master/R_data/train_set_imp2.rds')
train_set_imp <- train_set_imp %>% group_by(ctry) %>% mutate(ctry_lpi = mean(lpi, na.rm = TRUE)) %>% ungroup()
train_set_imp <- train_set_imp %>% group_by(ctry) %>% mutate(ctry_Interpersonal_Index = mean(Interpersonal_Index, na.rm = TRUE)) %>% ungroup()
train_set_imp <- train_set_imp %>% group_by(ctry) %>% mutate(ctry_Institutional_Index = mean(Institutional_Index , na.rm = TRUE)) %>% ungroup()
country_means_reference1 <- train_set_imp %>% group_by(ctry) %>%
summarise(
ctry_lpi = mean(ctry_lpi, na.rm = TRUE),
ctry_Interpersonal_Index = mean(ctry_Interpersonal_Index, na.rm = TRUE),
ctry_Institutional_Index = mean(ctry_Institutional_Index , na.rm = TRUE)
)
#saveRDS(country_means_reference1, "country_means_reference2.rds")train_set_imp <- train_set_imp %>% select(-ctry)
pp_model <- preProcess(train_set_imp, method = c("center", "scale"))
#saveRDS(train_columns, "train_columns2.rds")
train_set_scaled <- predict(pp_model, newdata = train_set_imp)train_set_scaled$hesitancy <- ifelse(train_set_scaled$hesitancy == "hesitant", 1, 0)
train_dum <- fastDummies::dummy_columns(train_set_scaled, remove_selected_columns = TRUE)
X_train <- train_dum %>% select(-hesitancy)
train_columns <- colnames(X_train)
y_train <- train_dum$hesitancy
y_train <- ifelse(y_train == 1, "hesitant", "willing") ctrl <- trainControl(method = "cv", number = 5, classProbs = TRUE, summaryFunction = twoClassSummary, savePredictions = "final") set.seed(42)
GLM = train(x = as.data.frame(X_train), y = y_train, method = 'glm', family = 'binomial', trControl = ctrl) Warning in train.default(x = as.data.frame(X_train), y = y_train, method =
"glm", : The metric "Accuracy" was not in the result set. ROC will be used
instead.GLM$results$ROC[1] 0.6294632RF = train(x = as.data.frame(X_train), y = y_train, method = "ranger", num.trees = 800,
importance = "permutation", verbose = FALSE, trControl = ctrl, metric = "ROC")
saveRDS(RF, '/Users/PaulMUTAMBA/Documents/Old Pc/KUL QASS/Master/R_data/RF2.rds')RF_adult <- readRDS('/Users/PaulMUTAMBA/Documents/Old Pc/KUL QASS/Master/R_data/RF2.rds')
max(RF_adult$results$ROC)[1] 0.6736716Building the Function
train_tuned_xgboost <- function(train_data, save_path = "final_xgboost_model.rds", seed = 42) {
# We prepare the target variable
train_data$hesitancy <- factor(train_data$hesitancy, levels = c("hesitant", "willing"))
# We Build the Recipe
xgb_recipe <- recipe(hesitancy ~ ., data = train_data) %>%
step_center(all_numeric_predictors()) %>%
step_scale(all_numeric_predictors()) %>%
step_dummy(all_nominal_predictors(), one_hot = TRUE)
# we Define the Spec
xgb_spec <- boost_tree(
trees = tune(), # Let the model find the best number of trees
tree_depth = tune(), # Let the model find the best depth
learn_rate = tune(), # Let the model find the best learning rate (eta)
min_n = tune(), # Controls overfitting at the node level
loss_reduction = tune(), # Controls tree pruning (gamma)
mtry = tune()
) %>%
set_engine("xgboost") %>%
set_mode("classification")
# We Bundle into a Workflow
xgb_workflow <- workflow() %>%
add_recipe(xgb_recipe) %>%
add_model(xgb_spec)
# We Create CV Folds
set.seed(seed)
cv_folds <- vfold_cv(train_data, v = 5)
# --------------------------------------------------------------------------
# THE TUNING MAGIC
# --------------------------------------------------------------------------
cat("Tuning hyperparameters (Training 50 models: 10 combos * 5 folds)...\n")
cat("This might take a minute or two!\n")
set.seed(seed)
tune_results <- tune_grid(
xgb_workflow,
resamples = cv_folds,
grid = 30, # Automatically tries 10 smart combinations
metrics = metric_set(roc_auc), # Optimize for best AUC
control = control_grid(save_pred = TRUE)
)
# We extract the absolute best parameters
best_params <- select_best(tune_results, metric = "roc_auc")
cat("\nThe Best Tuning Parameters Found:\n")
print(best_params)
oof_preds <- collect_predictions(tune_results, parameters = best_params)
cat("\nAUC Check:\n")
cv_auc_val <- roc_auc(oof_preds, truth = hesitancy, .pred_hesitant)
print(cv_auc_val)
# We Find the Youden Threshold using 'probably'
threshold_data <- oof_preds %>%
threshold_perf(truth = hesitancy, estimate = .pred_hesitant, metrics = metric_set(j_index))
optimal_threshold <- threshold_data %>%
filter(.metric == "j_index") %>%
arrange(desc(.estimate)) %>%
slice(1) %>%
pull(.threshold)
cat("\nOptimal CV Threshold (Youden):", optimal_threshold, "\n")
cat("\nTop 5 models by AUC:\n")
print(show_best(tune_results, metric = "roc_auc"))
# --------------------------------------------------------------------------
# Finalize and Train
# --------------------------------------------------------------------------
cat("\nTraining the final model using the best parameters...\n")
# This automatically replaces the "tune()" placeholders with the winning numbers
final_workflow <- finalize_workflow(xgb_workflow, best_params)
# This fit the final model to the entire training set
final_fit <- fit(final_workflow, data = train_data)
# --------------------------------------------------------------------------
# Here we save the Model
# --------------------------------------------------------------------------
cat(paste0("\nSaving the final model to disk as: ", save_path, "\n"))
saveRDS(final_fit, file = save_path)
cat("Done! Your optimized model is ready.\n")
# Here we return the Results
return(list(
model = final_fit,
best_params = best_params,
optimal_threshold = optimal_threshold,
cv_auc = cv_auc_val$.estimate
))
}
##____FIN_____##Training the XGBoost Model
my_results <- train_tuned_xgboost(
train_data = train_set_imp,
save_path = "my_optimized_xgboost3.rds"
)xgbmodel <- readRDS('/Users/PaulMUTAMBA/Documents/Old Pc/KUL QASS/Master/Codes/my_optimized_xgboost3.rds')| Models | AUC |
|---|---|
| GLM | 0.6294 |
| RF | 0.6736 |
| XGBOOST | 0.686 |
train_set_imp_young <- train_set_imp %>% filter(age %in% c(18:35))
pp_model <- preProcess(train_set_imp_young, method = c("center", "scale"))
train_set_scaled <- predict(pp_model, newdata = train_set_imp_young)
train_set_scaled$hesitancy <- ifelse(train_set_scaled$hesitancy == "hesitant", 1, 0)
train_dum1 <- fastDummies::dummy_columns(train_set_scaled, remove_selected_columns = TRUE)
X_train <- train_dum1 %>% select(-hesitancy)
y_train <- train_dum1$hesitancy
y_train <- ifelse(y_train == 1, "hesitant", "willing")
set.seed(42)
GLM_young = train(x = as.data.frame(X_train), y = y_train, method = 'glm', family = 'binomial', trControl = ctrl)
saveRDS(GLM_young, '/Users/PaulMUTAMBA/Documents/Old Pc/KUL QASS/Master/R_data/GLM_young2.rds')GLM_young <- readRDS('/Users/PaulMUTAMBA/Documents/Old Pc/KUL QASS/Master/R_data/GLM_young2.rds')
max(GLM_young$results$ROC)[1] 0.6242315train_set_imp_young <- train_set_imp %>% filter(age %in% c(18:35))
pp_model <- preProcess(train_set_imp_young, method = c("center", "scale"))
train_set_scaled <- predict(pp_model, newdata = train_set_imp_young)
train_set_scaled$hesitancy <- ifelse(train_set_scaled$hesitancy == "hesitant", 1, 0)
train_dum1 <- fastDummies::dummy_columns(train_set_scaled, remove_selected_columns = TRUE)
X_train <- train_dum1 %>% select(-hesitancy)
y_train <- train_dum1$hesitancy
y_train <- ifelse(y_train == 1, "hesitant", "willing")
set.seed(42)
RF_young = train(x = as.data.frame(X_train), y = y_train , method = "ranger", num.trees = 800,
importance = "permutation", verbose = FALSE, trControl = ctrl, metric = "ROC")
saveRDS(RF_young, '/Users/PaulMUTAMBA/Documents/Old Pc/KUL QASS/Master/R_data/RF_young1.rds')RF_young <- readRDS('/Users/PaulMUTAMBA/Documents/Old Pc/KUL QASS/Master/R_data/RF_young1.rds')
max(RF_young$results$ROC)[1] 0.6613971train_set_imp_young <- train_set_imp %>% filter(age %in% c(18:35))
young_model <- train_tuned_xgboost(
train_data = train_set_imp_young,
save_path = "my_optimized_xgboost_young_subsample2.rds"
)| Models | AUC |
|---|---|
| GLM | 0.6242 |
| RF | 0.661 |
| XGBOOST | 0.674 |
train_set_imp_adult <- train_set_imp %>% filter(age > 35)
pp_model_adult <- preProcess(train_set_imp_adult, method = c("center", "scale"))
#saveRDS(pp_model_adult, "pp_model_adult.rds" )
train_set_scaled <- predict(pp_model_adult, newdata = train_set_imp_adult)
train_set_scaled$hesitancy <- ifelse(train_set_scaled$hesitancy == "hesitant", 1, 0)
train_dum1 <- fastDummies::dummy_columns(train_set_scaled, remove_selected_columns = TRUE)
X_train <- train_dum1 %>% select(-hesitancy)
y_train <- train_dum1$hesitancy
y_train <- ifelse(y_train == 1, "hesitant", "willing")
set.seed(42)
GLM_adult = train(x = as.data.frame(X_train), y = y_train, method = 'glm',
family = 'binomial', trControl = ctrl) Warning in train.default(x = as.data.frame(X_train), y = y_train, method =
"glm", : The metric "Accuracy" was not in the result set. ROC will be used
instead.GLM_adult$results$ROC[1] 0.6318756#saveRDS(GLM_adult, '/Users/PaulMUTAMBA/Documents/Old Pc/KUL QASS/Master/Codes/GLM_adult1.rds')train_set_imp_adult <- train_set_imp %>% filter(age > 35)
pp_model <- preProcess(train_set_imp_adult, method = c("center", "scale"))
train_set_scaled <- predict(pp_model, newdata = train_set_imp_adult)
train_set_scaled$hesitancy <- ifelse(train_set_scaled$hesitancy == "hesitant", 1, 0)
train_dum1 <- fastDummies::dummy_columns(train_set_scaled, remove_selected_columns = TRUE)
X_train <- train_dum1 %>% select(-hesitancy)
y_train <- train_dum1$hesitancy
y_train <- ifelse(y_train == 1, "hesitant", "willing")
set.seed(42)
RF_adult = train(x = as.data.frame(X_train), y = y_train , method = "ranger", num.trees = 800,
importance = "permutation", verbose = FALSE, trControl = ctrl, metric = "ROC")
saveRDS(RF_adult, '/Users/PaulMUTAMBA/Documents/Old Pc/KUL QASS/Master/R_data/RF_adult1.rds')RF_adult <- readRDS('/Users/PaulMUTAMBA/Documents/Old Pc/KUL QASS/Master/R_data/RF_adult1.rds')
max(RF_adult$results$ROC)[1] 0.6768533train_set_imp_adult <- train_set_imp %>% filter(age > 35)
Adult_model <- train_tuned_xgboost(
train_data = train_set_imp_adult,
save_path = "my_optimized_xgboost_adult_subsample2.rds"
)| Models | AUC |
|---|---|
| GLM | 0.6318 |
| RF | 0.676 |
| XGBOOST | 0.691 |
glm_oof_preds <- GLM_adult$pred
# 2. Calculate the ROC and Youden Index
# (Assuming 'hesitant' is the column containing the probability of hesitancy)
glm_roc_obj <- roc(
response = glm_oof_preds$obs,
predictor = glm_oof_preds$hesitant,
quiet = TRUE
)
glm_best_coords <- pROC::coords(
glm_roc_obj, x = "best", best.method = "youden",
ret = c("threshold", "specificity", "sensitivity")
)
# 3. Save the optimal threshold
glm_optimal_threshold <- glm_best_coords$threshold[1]
cat("GLM Optimal CV Threshold (Youden):", glm_optimal_threshold, "\n")GLM Optimal CV Threshold (Youden): 0.5712388 cat("CV Specificity:", glm_best_coords$specificity[1], "| CV Sensitivity:", glm_best_coords$sensitivity[1], "\n")CV Specificity: 0.5796634 | CV Sensitivity: 0.6146026 Here we repeat the preprocess steps applied for the train_set.
test_set <- readRDS('/Users/PaulMUTAMBA/Documents/Old Pc/KUL QASS/Master/R_data/test_set.rds')
test_set <- test_set %>% mutate(across(where(~class(.) == "character"), ~as.factor(.)))
diff_cols <- c("difficulty_schl_services", "difficulty_hsptl_services",
"difficulty_id_services", "difficulty_police_services")
respect_cols <- c("feel_rspect_schl", "feel_rspect_hsptl", "feel_rspect_id")
df <- test_set %>%
mutate(
# Overall Service Utilization (1 = Yes, 0 = No)
# Checks if any of the columns are NOT "no contact/answer"
service_utilization = as.integer(if_any(all_of(diff_cols), ~ . != "no contact/answer" & !is.na(.))),
# Convert text to numbers for Difficulty (creates temporary columns)
across(all_of(diff_cols),
~ case_when(
. == "very easy" ~ 1,
. == "easy" ~ 2,
. == "difficult" ~ 3,
. == "very difficult" ~ 4,
TRUE ~ NA_real_ # Everything else, including "no contact", becomes NA
), .names = "{.col}_num"),
# Convert text to numbers for Respect (creates temporary columns)
across(all_of(respect_cols),
~ case_when(
. == "not at all" ~ 1,
. == "a little bit" ~ 2,
. == "somewhat" ~ 3,
. == "a lot" ~ 4,
TRUE ~ NA_real_
), .names = "{.col}_num")
) %>%
# Calculate the row means, ignoring NAs
rowwise() %>%
mutate(
avg_difficulty = mean(c_across(ends_with("_num") & contains("difficulty")), na.rm = TRUE),
avg_respect = mean(c_across(ends_with("_num") & contains("rspect")), na.rm = TRUE)
) %>%
ungroup() %>%
# Drop the temporary numeric columns to keep the dataframe clean
select(-ends_with("_num"), -starts_with("diff"), -starts_with("feel"))
df <- df %>%
mutate(
avg_diff_model = if_else(is.na(avg_difficulty), 0, avg_difficulty),
avg_resp_model = if_else(is.na(avg_respect), 0, avg_respect)
) %>% select(-avg_difficulty, -avg_respect)
df_mlvl <- df
df_mlvl$service_utilization <- ifelse(df_mlvl$service_utilization == "yes", 1, 0)
col_num <- colnames(df_mlvl %>% select_if(~class(.) == "numeric"))
df_mlvl$Interpersonal_Index <- rowMeans(df_mlvl[, c("trst_relative", "trst_neighbours", "trst_people_other_religion", "trst_other_people_know")], na.rm=TRUE)
df_mlvl$Institutional_Index <- rowMeans(df_mlvl[, c("trst_pres", "trst_police", "trst_courts")], na.rm=TRUE)
df_mlvl <- df_mlvl %>%
mutate(
language_profile = case_when(
lang_home == "colonial" & lang_interview == "colonial" ~ "Consistent Colonial",
lang_home == "indigeneous/others" & lang_interview == "indigeneous/others" ~ "Consistent Indigenous",
(lang_home == "colonial" & lang_interview == "indigeneous/others") |
(lang_home == "indigeneous/others" & lang_interview == "colonial") ~ "Mixed/Code-Switcher",
TRUE ~ NA_character_
)
)
civic_cols <- c("attnd_cmty_meeting", "raise_issue", "protest_march")
df_mlvl <- df_mlvl %>%
mutate(
# Convert text to numbers using case_when (creates temporary columns)
across(all_of(civic_cols),
~ case_when(
. == "Active" ~ 3,
. == "Willing" ~ 2,
. == "Resistant" ~ 1,
TRUE ~ NA_real_ # "No Answer" and anything else becomes NA
), .names = "{.col}_num")
) %>%
# Calculate the row means, ignoring NAs
rowwise() %>%
mutate(
civic_participation_score = mean(c_across(ends_with("_num") & contains(c("attnd", "raise", "protest"))), na.rm = TRUE)
) %>%
ungroup() %>%
# Clean up by dropping the temporary numeric columns
select(-ends_with("_num"))
df_mlvl$lang_policy <- fct_collapse(
df_mlvl$indigenous_lang_integration,
"no integration" = "no integration",
"limited integration" = "limited integration",
"moderate or extensive" = c("moderate integration", "extensive integration")
)
test_set <- df_mlvl %>% mutate(across(where(~class(.) == "character"), ~as.factor(.))) %>%
select(hesitancy, language_profile, avg_diff_model, avg_resp_model, service_utilization,
Interpersonal_Index, Institutional_Index, civic_participation_score, gndr, age, urban_rural, lvl_edu, religion,
lpi,ctry, lang_policy)Next we applied similar imputation process to the test_set.
##----imputation------##
unique_countries <- unique(test_set$ctry)
test_set$ctry <- as.factor(test_set$ctry)
country_dfs <- list()
for (country in unique_countries) {
country_df <- test_set[test_set$ctry == country, ]
country_dfs[[country]] <- as.data.frame(country_df)
}
library(mi)
library(doParallel)
library(missForest)
#summary(train_set$avg_diff_model)
fact_var <- colnames(test_set %>% select_if(~class(.) %in% c("character", "factor")))
num_var <- colnames(test_set %>% select_if(~!class(.) %in% c("character", "factor")))
mdf <- list()
for (country in unique_countries) {
mdf[[country]] <- missing_data.frame (as.data.frame(country_dfs[[country]]) )
mdf[[country]] <- change (mdf[[country]] , y = num_var, what = "type", to = "continuous")
mdf[[country]] <- change (mdf[[country]] , y = fact_var, what = "type", to = "ordered-categorical")
}
##-----MissForest Imputation----##
registerDoParallel(cores=4)
set.seed(42)
imputed_missForests <- list()
for (country in unique_countries){
imputed_missForests[[country]] <- missForest(country_dfs[[country]], verbose = TRUE, parallelize = "forests")
}
imputed_missForests[[1]]$OOBerror[1]
missForest_imputed <- data.frame(
original = country_dfs[[1]]$trst_courts,
imputed_missForest = imputed_missForests[[1]]$ximp$trst_courts
)
test_set_imp <- data.frame()
for (country in unique_countries){
test_set_imp <- rbind(test_set_imp, as.data.frame(imputed_missForests[[country]]$ximp))
}
test_set_imp <- test_set_imp %>% mutate(across(where(is.numeric), round))
#saveRDS(test_set_imp, "test_set_imp.rds")Country Level Variables [From train_set]
Given that the national-level indices were calculated using the training datasets, which contain more variables, it is recommended to use the national-level indices derived from the training dataset rather than creating new ones from the test dataset. This will not affect data leakage as such, as all instances from the same country have the same national-level scores, including those that the training model has yet to encounter.
test_set_imp <- readRDS('/Users/PaulMUTAMBA/Documents/Old Pc/KUL QASS/Master/Codes/test_set_imp.rds')
##----Country Level Variable for test-----##
country_means_reference1 <- readRDS("/Users/PaulMUTAMBA/Test/country_means_reference2.rds")
test_set_imp <- test_set_imp %>%
left_join(country_means_reference1, by = "ctry")
test_set_imp <- test_set_imp %>% select(-ctry)model_full <- readRDS('/Users/PaulMUTAMBA/Documents/Old Pc/KUL QASS/Master/Codes/my_optimized_xgboost3.rds')
test_results <- workflows:::predict.workflow(model_full, new_data = test_set_imp, type = "prob") %>%
bind_cols(test_set_imp %>% select(hesitancy))
test_results <- test_results %>%
mutate(
final_prediction = factor(
ifelse(.pred_hesitant >= 0.55, "hesitant", "willing"),
levels = c("hesitant", "willing")
)
)test_results$hesitancy <- factor(test_results$hesitancy, levels = c("hesitant", "willing"))
cat("\n--- Final Test Set Confusion Matrix ---\n")
--- Final Test Set Confusion Matrix ---confusionMatrix(
data = test_results$final_prediction,
reference = test_results$hesitancy,
positive = "hesitant"
)Confusion Matrix and Statistics
Reference
Prediction hesitant willing
hesitant 1725 793
willing 822 1180
Accuracy : 0.6427
95% CI : (0.6285, 0.6567)
No Information Rate : 0.5635
P-Value [Acc > NIR] : <2e-16
Kappa : 0.2749
Mcnemar's Test P-Value : 0.486
Sensitivity : 0.6773
Specificity : 0.5981
Pos Pred Value : 0.6851
Neg Pred Value : 0.5894
Prevalence : 0.5635
Detection Rate : 0.3816
Detection Prevalence : 0.5571
Balanced Accuracy : 0.6377
'Positive' Class : hesitant
library("yardstick")
cat("\n--- Final Test Set AUC ---\n")
--- Final Test Set AUC ---test_auc <- test_results %>%
roc_auc(truth = hesitancy, .pred_hesitant)
print(test_auc)# A tibble: 1 × 3
.metric .estimator .estimate
<chr> <chr> <dbl>
1 roc_auc binary 0.687SHAP
xgb_raw_model <- extract_fit_engine(model_full)
test_recipe <- extract_recipe(model_full)
baked_test <- bake(test_recipe, new_data = test_set_imp)
X_test_matrix <- baked_test %>%
select(-hesitancy) %>%
as.matrix()
shp <- shapviz(xgb_raw_model, X_pred = X_test_matrix)
p1 <- sv_importance(shp, kind = "beeswarm")
print(p1)
p2 <- sv_importance(shp, kind = "bar")
print(p2)
Permutation Importance
set.seed(42)
perm_importance_plot <- model_full %>%
extract_fit_parsnip() %>%
vip(
method = "permute",
train = baked_test,
target = "hesitancy",
metric = "roc_auc",
reference_class = "hesitant",
event_level = "first",
pred_wrapper = function(object, newdata) {
# 1. Generate the predictions
preds <- predict(object, newdata, type = "prob")
# 2. Check the format and return accordingly
if (is.data.frame(preds)) {
# If tidymodels returns a dataframe, grab the specific column
return(preds$.pred_hesitant)
} else {
# If the raw engine returns a numeric vector, just return it
return(as.numeric(preds))
}
},
nsim = 5,
geom = "boxplot",
num_features = 20,
aesthetics = list(fill = "darkorange")
) +
theme_minimal() +
labs(
title = NULL,
subtitle = NULL,
x = "Feature",
y = "Decrease in AUC"
)
print(perm_importance_plot)
test_set_imp_yg <- test_set_imp %>% filter(age %in% c(18:35))
model_young <- readRDS('/Users/PaulMUTAMBA/Documents/Old Pc/KUL QASS/Master/Codes/my_optimized_xgboost_young_subsample2.rds')
test_results <- workflows:::predict.workflow(model_young, new_data = test_set_imp_yg, type = "prob") %>%
bind_cols(test_set_imp_yg %>% select(hesitancy))
# 3. Apply your optimal Youden threshold
# (Assuming your threshold from the previous script is saved in 'optimal_threshold')
test_results <- test_results %>%
mutate(
# If the probability of being hesitant is higher than the threshold, classify as hesitant
final_prediction = factor(
ifelse(.pred_hesitant >= 0.55, "hesitant", "willing"),
levels = c("hesitant", "willing")
)
)test_results$hesitancy <- factor(test_results$hesitancy, levels = c("hesitant", "willing"))
# --------------------------------------------------------------------------
# 4. Evaluate the Results
# --------------------------------------------------------------------------
cat("\n--- Final Test Set Confusion Matrix ---\n")
--- Final Test Set Confusion Matrix ---confusionMatrix(
data = test_results$final_prediction,
reference = test_results$hesitancy,
positive = "hesitant"
)Confusion Matrix and Statistics
Reference
Prediction hesitant willing
hesitant 943 457
willing 516 680
Accuracy : 0.6252
95% CI : (0.6062, 0.6439)
No Information Rate : 0.562
P-Value [Acc > NIR] : 3.744e-11
Kappa : 0.243
Mcnemar's Test P-Value : 0.06297
Sensitivity : 0.6463
Specificity : 0.5981
Pos Pred Value : 0.6736
Neg Pred Value : 0.5686
Prevalence : 0.5620
Detection Rate : 0.3633
Detection Prevalence : 0.5393
Balanced Accuracy : 0.6222
'Positive' Class : hesitant
cat("\n--- Final Test Set AUC ---\n")
--- Final Test Set AUC ---test_auc <- test_results %>%
roc_auc(truth = hesitancy, .pred_hesitant)
print(test_auc)# A tibble: 1 × 3
.metric .estimator .estimate
<chr> <chr> <dbl>
1 roc_auc binary 0.674SHAP
xgb_raw_model <- extract_fit_engine(model_young)
test_recipe <- extract_recipe(model_young)
baked_test <- bake(test_recipe, new_data = test_set_imp_yg)
X_test_matrix <- baked_test %>%
select(-hesitancy) %>%
as.matrix()
shp <- shapviz(xgb_raw_model, X_pred = X_test_matrix)
p1 <- sv_importance(shp, kind = "beeswarm")
print(p1)
p2 <- sv_importance(shp, kind = "bar")
print(p2)
Permutation Importance
set.seed(42)
perm_importance_plot <- model_young %>%
extract_fit_parsnip() %>% # Extract the model engine
vip(
method = "permute",
train = baked_test, # Using the baked data (67 columns)
target = "hesitancy",
metric = "roc_auc",
reference_class = "hesitant",
event_level = "first",
pred_wrapper = function(object, newdata) {
# 1. Generate the predictions
preds <- predict(object, newdata, type = "prob")
# 2. Check the format and return accordingly
if (is.data.frame(preds)) {
# If tidymodels returns a dataframe, grab the specific column
return(preds$.pred_hesitant)
} else {
# If the raw engine returns a numeric vector, just return it
return(as.numeric(preds))
}
},
nsim = 5,
geom = "boxplot",
num_features = 20,
aesthetics = list(fill = "darkblue")
) +
theme_minimal() +
labs(
title = "Permutation Importance (Youths Between 18 and 35)",
subtitle = "Drop in ROC AUC when a feature's values are randomly shuffled",
x = "Decrease in AUC",
y = "Feature"
)
print(perm_importance_plot)
test_set_imp_adult <- test_set_imp %>% filter(age > 35)
model_adult <- readRDS('/Users/PaulMUTAMBA/Documents/Old Pc/KUL QASS/Master/Codes/my_optimized_xgboost_adult_subsample2.rds')
test_results <- workflows:::predict.workflow(model_adult, new_data = test_set_imp_adult, type = "prob") %>%
bind_cols(test_set_imp_adult %>% select(hesitancy))
test_results <- test_results %>%
mutate(
final_prediction = factor(
ifelse(.pred_hesitant >= 0.575, "hesitant", "willing"),
levels = c("hesitant", "willing")
)
)test_results$hesitancy <- factor(test_results$hesitancy, levels = c("hesitant", "willing"))
# --------------------------------------------------------------------------
# 4. Evaluate the Results
# --------------------------------------------------------------------------
cat("\n--- Final Test Set Confusion Matrix ---\n")
--- Final Test Set Confusion Matrix ---confusionMatrix(
data = test_results$final_prediction,
reference = test_results$hesitancy,
positive = "hesitant"
)Confusion Matrix and Statistics
Reference
Prediction hesitant willing
hesitant 694 294
willing 394 542
Accuracy : 0.6424
95% CI : (0.6205, 0.6639)
No Information Rate : 0.5655
P-Value [Acc > NIR] : 3.856e-12
Kappa : 0.2823
Mcnemar's Test P-Value : 0.0001604
Sensitivity : 0.6379
Specificity : 0.6483
Pos Pred Value : 0.7024
Neg Pred Value : 0.5791
Prevalence : 0.5655
Detection Rate : 0.3607
Detection Prevalence : 0.5135
Balanced Accuracy : 0.6431
'Positive' Class : hesitant
cat("\n--- Final Test Set AUC ---\n")
--- Final Test Set AUC ---test_auc <- test_results %>%
roc_auc(truth = hesitancy, .pred_hesitant)
print(test_auc)# A tibble: 1 × 3
.metric .estimator .estimate
<chr> <chr> <dbl>
1 roc_auc binary 0.697SHAP
xgb_raw_model <- extract_fit_engine(model_adult)test_recipe <- extract_recipe(model_adult)
baked_test <- bake(test_recipe, new_data = test_set_imp_adult)
X_test_matrix <- baked_test %>%
select(-hesitancy) %>%
as.matrix()
shp <- shapviz(xgb_raw_model, X_pred = X_test_matrix)
p1 <- sv_importance(shp, kind = "beeswarm")
print(p1)
p2 <- sv_importance(shp, kind = "bar")
print(p2)
Permutation
set.seed(42)
perm_importance_plot <- model_adult %>%
extract_fit_parsnip() %>% # Extract the model engine
vip(
method = "permute",
train = baked_test, # Using the baked data (67 columns)
target = "hesitancy",
metric = "roc_auc",
reference_class = "hesitant",
event_level = "first",
pred_wrapper = function(object, newdata) {
# 1. Generate the predictions
preds <- predict(object, newdata, type = "prob")
# 2. Check the format and return accordingly
if (is.data.frame(preds)) {
# If tidymodels returns a dataframe, grab the specific column
return(preds$.pred_hesitant)
} else {
# If the raw engine returns a numeric vector, just return it
return(as.numeric(preds))
}
},
nsim = 5,
geom = "boxplot",
num_features = 20,
aesthetics = list(fill = "darkorange")
) +
theme_minimal() +
labs(
title = "Permutation Importance (Adult > 35)",
subtitle = "Drop in ROC AUC when a feature's values are randomly shuffled",
x = "Decrease in AUC",
y = "Feature"
)
print(perm_importance_plot)
Similarly to the Machine Learning Process, we have employed advanced data processing techniques here. We have converted categorical variables into numerical variables and combined various scores into a single score. This mainly applies to indices relating to institutions and interpersonal relationships, as well as scores for civic engagement. Here we call the df_ml dataset saved before the machine learning steps.
df <- readRDS("/Users/PaulMUTAMBA/Documents/Old Pc/KUL QASS/Master/R_data/df_ml.rds")
#---- Merge Difficulties and Respect-------##
diff_cols <- c("difficulty_schl_services", "difficulty_hsptl_services",
"difficulty_id_services", "difficulty_police_services")
respect_cols <- c("feel_rspect_schl", "feel_rspect_hsptl", "feel_rspect_id")
df <- df %>%
mutate(
service_utilization = as.integer(if_any(all_of(diff_cols), ~ . != "no contact/answer" & !is.na(.))),
across(all_of(diff_cols),
~ case_when(
. == "very easy" ~ 1,
. == "easy" ~ 2,
. == "difficult" ~ 3,
. == "very difficult" ~ 4,
TRUE ~ NA_real_
), .names = "{.col}_num"),
across(all_of(respect_cols),
~ case_when(
. == "not at all" ~ 1,
. == "a little bit" ~ 2,
. == "somewhat" ~ 3,
. == "a lot" ~ 4,
TRUE ~ NA_real_
), .names = "{.col}_num")
) %>%
rowwise() %>%
mutate(
avg_difficulty = mean(c_across(ends_with("_num") & contains("difficulty")), na.rm = TRUE),
avg_respect = mean(c_across(ends_with("_num") & contains("rspect")), na.rm = TRUE)
) %>%
ungroup() %>%
select(-ends_with("_num"), -starts_with("diff"), -starts_with("feel"))
df <- df %>%
mutate(
avg_diff_model = if_else(is.na(avg_difficulty), 0, avg_difficulty),
avg_resp_model = if_else(is.na(avg_respect), 0, avg_respect)
) %>% select(-avg_difficulty, -avg_respect)
df_ml <- df
##-----country personal trust------##
afro_perso <- df_ml %>%
mutate(personal_trust = rowMeans(pick(trst_relative, trst_neighbours, trst_people_other_religion, trst_other_people_know), na.rm = TRUE)) %>%
group_by(ctry) %>%
mutate(ctry_personal_trust = mean(personal_trust, na.rm = TRUE)) %>%
ungroup()
##------country lpi-------##
afro_lpi <- afro_perso %>% group_by(ctry) %>% mutate(ctry_lpi = mean(lpi, na.rm = TRUE)) %>% ungroup()
##------Country Levels Institutional Trust--------##
afro_lpi <- afro_lpi %>% mutate(mean_institutional_trus = rowMeans(pick(trst_pres, trst_police, trst_courts), na.rm = TRUE)) %>%
group_by(ctry) %>% mutate(country_intst_trst = mean(mean_institutional_trus, na.rm = TRUE)) %>% ungroup()
##---------Multi level analysis---------###
##------To augment interation to have a model that converges--------##
afro_lpi$service_utilization <- as.factor(ifelse(afro_lpi$service_utilization == 1, "yes", "no"))
df_mlvl <- afro_lpi %>% ungroup()
#df_mlvl <- df_mlvl %>% mutate(across(where(is.numeric), ~as.numeric(scale(.))))
df_mlvl$service_utilization <- ifelse(df_mlvl$service_utilization == "yes", 1, 0)
col_num <- colnames(df_mlvl %>% select_if(~class(.) == "numeric"))
df_mlvl <- df_mlvl %>% select(-c("gdp", "gini"))
df_mlvl <- df_mlvl %>%
mutate(across(where(is.character), as.factor))
df_mlvl <- df_mlvl %>%
mutate(
across(
c(attnd_cmty_meeting, raise_issue, protest_march),
~ .x %>%
tolower() %>%
na_if("no answer") %>%
factor()
)
)Similarly to the ML Process, we apply the same imputation techniques for the Multilevel Analysis.
unique_countries <- unique(df_mlvl$ctry)
df_mlvl$ctry <- as.factor(df_mlvl$ctry)
country_dfs <- list()
for (country in unique_countries) {
country_df <- df_mlvl[df_mlvl$ctry == country, ]
country_dfs[[country]] <- as.data.frame(country_df)
}
fact_var <- colnames(df_mlvl %>% select_if(~class(.) %in% c("character", "factor")))
num_var <- colnames(df_mlvl %>% select_if(~!class(.) %in% c("character", "factor")))
mdf <- list()
for (country in unique_countries) {
mdf[[country]] <- missing_data.frame (as.data.frame(country_dfs[[country]]) )
mdf[[country]] <- change (mdf[[country]] , y = num_var, what = "type", to = "continuous")
mdf[[country]] <- change (mdf[[country]] , y = fact_var, what = "type", to = "ordered-categorical")
}
registerDoParallel(cores=4)
set.seed(42)
imputed_missForests <- list()
for (country in unique_countries){
imputed_missForests[[country]] <- missForest(country_dfs[[country]], verbose = TRUE, parallelize = "forests")
}
imputed_missForests[[1]]$OOBerror[1]
df_mlvl_imp <- data.frame()
for (country in unique_countries){
df_mlvl_imp <- rbind(df_mlvl_imp, as.data.frame(imputed_missForests[[country]]$ximp))
}
df_mlvl_imp <- df_mlvl_imp %>% mutate(across(where(is.numeric), round))
sum(is.na(df_mlvl_imp))
saveRDS(df_mlvl_imp, "/Users/PaulMUTAMBA/Documents/Old Pc/KUL QASS/Master/R_data/df_mlvl_imp.rds")Although this procedure could have been applied during the machine learning stage, we felt it would be redundant, as its sole purpose is to provide insight into the statistical significance of aggregating institutional and interpersonal scores. That is why we present it here, during the multilevel analysis, where this aggregation is of particular importance.
df_mlvl_imp <- readRDS("/Users/PaulMUTAMBA/Documents/Old Pc/KUL QASS/Master/R_data/df_mlvl_imp.rds")
trust_df <- df_mlvl_imp %>% select(starts_with("trst"))
fa.parallel(trust_df, fa="both", n.iter=100) 
Parallel analysis suggests that the number of factors = 3 and the number of components = 2 fa.diagram(fa(trust_df, nfactors=3, rotate="none", fm="ml"))
pca <- prcomp(trust_df, center = TRUE, scale. = TRUE)
fviz_pca_var(pca)
fviz_eig(pca, addlabels = TRUE, ylim = c(0, 60))
pca$rotation PC1 PC2 PC3 PC4
trst_pres 0.3292325 -0.4049075 0.14928384 0.83650448
trst_police 0.3562629 -0.4671808 -0.14040591 -0.33201586
trst_courts 0.3230553 -0.5261637 0.07390202 -0.39095892
trst_relative 0.3688159 0.2566377 -0.71641825 0.15061510
trst_neighbours 0.4417233 0.2615110 -0.23281502 -0.07644107
trst_people_other_religion 0.3825865 0.3352032 0.51878269 -0.04952017
trst_other_people_know 0.4275522 0.3071346 0.34050784 -0.07871758
PC5 PC6 PC7
trst_pres -0.06356382 -0.01176705 0.03810039
trst_police 0.29084765 -0.65894506 0.07790748
trst_courts -0.13992830 0.65286624 -0.12108821
trst_relative 0.36083853 0.21923576 -0.28966040
trst_neighbours -0.56736100 -0.03863639 0.59418447
trst_people_other_religion 0.57611103 0.17649737 0.32653237
trst_other_people_know -0.32829938 -0.24230190 -0.65896488pca$sdev[1] 1.7200766 1.1836108 0.8573255 0.7845719 0.6872963 0.6686986 0.6085292alpha1test <- psych::alpha(df_mlvl_imp[,c("trst_relative", "trst_neighbours", "trst_people_other_religion", "trst_other_people_know")])Warning in response.frequencies(x, max = max): response.frequency has been
deprecated and replaced with responseFrequecy. Please fix your callalpha1test$alpha.drop$raw_alpha[1] 0.7700963 0.6969408 0.7443662 0.7006797alpha1test$item.stats$r.drop[1] 0.5035220 0.6491901 0.5611146 0.6442620alpha1test$item.stats$r.cor[1] 0.5748357 0.7388823 0.6407304 0.7315999df_mlvl_imp$Interpersonal_Index <- rowMeans(df_mlvl_imp[, c("trst_relative", "trst_neighbours", "trst_people_other_religion", "trst_other_people_know")], na.rm=TRUE)alpha2test<- psych::alpha(df_mlvl_imp[,c("trst_pres", "trst_police", "trst_courts")])Warning in response.frequencies(x, max = max): response.frequency has been
deprecated and replaced with responseFrequecy. Please fix your calldf_mlvl_imp$Institutional_Index <- rowMeans(df_mlvl_imp[, c("trst_pres", "trst_police", "trst_courts")], na.rm=TRUE)
alpha2test$alpha.drop$raw_alpha[1] 0.7026970 0.5898673 0.5929705alpha2test$item.stats$r.cor[1] 0.5581727 0.6801184 0.6779330df_mlvl <- df_mlvl_imp %>%
mutate(
language_profile = case_when(
lang_home == "colonial" & lang_interview == "colonial" ~ "Consistent Colonial",
lang_home == "indigeneous/others" & lang_interview == "indigeneous/others" ~ "Consistent Indigenous",
(lang_home == "colonial" & lang_interview == "indigeneous/others") |
(lang_home == "indigeneous/others" & lang_interview == "colonial") ~ "Mixed/Code-Switcher",
TRUE ~ NA_character_
)
)civic_cols <- c("attnd_cmty_meeting", "raise_issue", "protest_march")
df_mlvl <- df_mlvl %>%
mutate(
# 1. Convert text to numbers using case_when (creates temporary columns)
across(all_of(civic_cols),
~ case_when(
. == "active" ~ 3,
. == "willing" ~ 2,
TRUE ~ 1
), .names = "{.col}_num")
) %>%
# 2. Calculate the row means, ignoring NAs
rowwise() %>%
mutate(
civic_participation_score = mean(c_across(ends_with("_num") & contains(c("attnd", "raise", "protest"))), na.rm = TRUE)
) %>%
ungroup() %>%
# 3. Clean up by dropping the temporary numeric columns
select(-ends_with("_num"))df_mlvl <- df_mlvl %>% group_by(ctry) %>% mutate(ctry_lpi = mean(lpi, na.rm = TRUE)) %>% ungroup()
df_mlvl <- df_mlvl %>% group_by(ctry) %>% mutate(ctry_Interpersonal_Index = mean(Interpersonal_Index, na.rm = TRUE)) %>% ungroup()
df_mlvl <- df_mlvl %>% group_by(ctry) %>% mutate(ctry_Institutional_Index = mean(Institutional_Index , na.rm = TRUE)) %>% ungroup()
##____Mundlak Expression____##
df_mlvl <- df_mlvl %>%
mutate(
Interpersonal_Index_dm = Interpersonal_Index - ctry_Interpersonal_Index,
Institutional_Index_dm = Institutional_Index - ctry_Institutional_Index,
lpi_dm = lpi - ctry_lpi
)
saveRDS(df_mlvl, "/Users/PaulMUTAMBA/Documents/Old Pc/KUL QASS/Master/R_data/df_mlvl.rds")Following a thorough data cleansing process, we examine here the distribution of the various covariates across the outcome variables. First we confirm that we don’t have NA’s in our dataset.
df_mlvl <- readRDS("/Users/PaulMUTAMBA/Documents/Old Pc/KUL QASS/Master/R_data/df_mlvl.rds")
print(sapply(df_mlvl, function(x) sum(is.na(x)))) lang_home lang_interview
0 0
trst_pres trst_police
0 0
trst_courts trst_relative
0 0
trst_neighbours trst_people_other_religion
0 0
trst_other_people_know attnd_cmty_meeting
0 0
raise_issue protest_march
0 0
gndr age
0 0
urban_rural lvl_edu
0 0
religion lpi
0 0
indigenous_lang_integration hesitancy
0 0
ctry service_utilization
0 0
avg_diff_model avg_resp_model
0 0
personal_trust ctry_personal_trust
0 0
ctry_lpi mean_institutional_trus
0 0
country_intst_trst Interpersonal_Index
0 0
Institutional_Index language_profile
0 0
civic_participation_score ctry_Interpersonal_Index
0 0
ctry_Institutional_Index Interpersonal_Index_dm
0 0
Institutional_Index_dm lpi_dm
0 0 table1::table1(~. |hesitancy, df_mlvl[, c("hesitancy","language_profile", "service_utilization", "avg_diff_model",
"avg_resp_model", "Interpersonal_Index", "Institutional_Index",
"civic_participation_score", "indigenous_lang_integration", "ctry_Interpersonal_Index",
"ctry_Institutional_Index", "gndr", "age", "urban_rural",
"lvl_edu", "religion", "lpi", "ctry_lpi")])Registered S3 method overwritten by 'table1':
method from
[.indexed plyr| willing (N=9820) |
hesitant (N=12778) |
Overall (N=22598) |
|
|---|---|---|---|
| language_profile | |||
| Consistent Colonial | 1012 (10.3%) | 2000 (15.7%) | 3012 (13.3%) |
| Consistent Indigenous | 7044 (71.7%) | 7819 (61.2%) | 14863 (65.8%) |
| Mixed/Code-Switcher | 1764 (18.0%) | 2959 (23.2%) | 4723 (20.9%) |
| service_utilization | |||
| Mean (SD) | 0.727 (0.445) | 0.706 (0.455) | 0.715 (0.451) |
| Median [Min, Max] | 1.00 [0, 1.00] | 1.00 [0, 1.00] | 1.00 [0, 1.00] |
| avg_diff_model | |||
| Mean (SD) | 1.76 (1.31) | 1.75 (1.34) | 1.75 (1.33) |
| Median [Min, Max] | 2.00 [0, 4.00] | 2.00 [0, 4.00] | 2.00 [0, 4.00] |
| avg_resp_model | |||
| Mean (SD) | 2.09 (1.57) | 1.91 (1.54) | 1.99 (1.56) |
| Median [Min, Max] | 2.00 [0, 4.00] | 2.00 [0, 4.00] | 2.00 [0, 4.00] |
| Interpersonal_Index | |||
| Mean (SD) | 2.95 (0.693) | 2.81 (0.738) | 2.87 (0.722) |
| Median [Min, Max] | 3.00 [1.00, 4.00] | 3.00 [1.00, 4.00] | 3.00 [1.00, 4.00] |
| Institutional_Index | |||
| Mean (SD) | 2.58 (0.885) | 2.29 (0.877) | 2.41 (0.892) |
| Median [Min, Max] | 2.67 [1.00, 4.00] | 2.33 [1.00, 4.00] | 2.33 [1.00, 4.00] |
| civic_participation_score | |||
| Mean (SD) | 2.03 (0.489) | 1.95 (0.535) | 1.99 (0.517) |
| Median [Min, Max] | 2.00 [1.00, 3.00] | 2.00 [1.00, 3.00] | 2.00 [1.00, 3.00] |
| indigenous_lang_integration | |||
| no integration | 3346 (34.1%) | 4670 (36.5%) | 8016 (35.5%) |
| limited integration | 3638 (37.0%) | 4332 (33.9%) | 7970 (35.3%) |
| moderate integration | 1959 (19.9%) | 3034 (23.7%) | 4993 (22.1%) |
| extensive integration | 877 (8.9%) | 742 (5.8%) | 1619 (7.2%) |
| ctry_Interpersonal_Index | |||
| Mean (SD) | 2.89 (0.282) | 2.86 (0.306) | 2.87 (0.296) |
| Median [Min, Max] | 2.84 [1.97, 3.44] | 2.84 [1.97, 3.44] | 2.84 [1.97, 3.44] |
| ctry_Institutional_Index | |||
| Mean (SD) | 2.45 (0.391) | 2.38 (0.386) | 2.41 (0.390) |
| Median [Min, Max] | 2.41 [1.69, 3.31] | 2.40 [1.69, 3.31] | 2.41 [1.69, 3.31] |
| gndr | |||
| man | 4833 (49.2%) | 6357 (49.7%) | 11190 (49.5%) |
| woman | 4987 (50.8%) | 6421 (50.3%) | 11408 (50.5%) |
| age | |||
| Mean (SD) | 35.8 (14.3) | 36.0 (14.0) | 35.9 (14.1) |
| Median [Min, Max] | 32.0 [18.0, 100] | 33.0 [18.0, 112] | 33.0 [18.0, 112] |
| urban_rural | |||
| rural | 5921 (60.3%) | 6203 (48.5%) | 12124 (53.7%) |
| urban | 3899 (39.7%) | 6575 (51.5%) | 10474 (46.3%) |
| lvl_edu | |||
| None or Informal | 1946 (19.8%) | 2288 (17.9%) | 4234 (18.7%) |
| Primary | 3113 (31.7%) | 3455 (27.0%) | 6568 (29.1%) |
| Secondary | 3498 (35.6%) | 4723 (37.0%) | 8221 (36.4%) |
| Higher Education | 1246 (12.7%) | 2283 (17.9%) | 3529 (15.6%) |
| No Answer | 17 (0.2%) | 29 (0.2%) | 46 (0.2%) |
| religion | |||
| None/Traditional | 502 (5.1%) | 872 (6.8%) | 1374 (6.1%) |
| Muslim | 3289 (33.5%) | 3930 (30.8%) | 7219 (31.9%) |
| Protestant/Christian | 4215 (42.9%) | 5717 (44.7%) | 9932 (44.0%) |
| Other/No Answer | 256 (2.6%) | 380 (3.0%) | 636 (2.8%) |
| Catholic | 1558 (15.9%) | 1879 (14.7%) | 3437 (15.2%) |
| lpi | |||
| Mean (SD) | 2.51 (0.954) | 2.52 (0.989) | 2.52 (0.974) |
| Median [Min, Max] | 3.00 [1.00, 5.00] | 3.00 [1.00, 5.00] | 3.00 [1.00, 5.00] |
| ctry_lpi | |||
| Mean (SD) | 2.52 (0.302) | 2.52 (0.354) | 2.52 (0.332) |
| Median [Min, Max] | 2.46 [1.28, 3.23] | 2.49 [1.28, 3.23] | 2.49 [1.28, 3.23] |
table1::table1(~ctry|hesitancy, df_mlvl)| willing (N=9820) |
hesitant (N=12778) |
Overall (N=22598) |
|
|---|---|---|---|
| ctry | |||
| angola | 89 (0.9%) | 61 (0.5%) | 150 (0.7%) |
| benin | 286 (2.9%) | 289 (2.3%) | 575 (2.5%) |
| botswana | 48 (0.5%) | 46 (0.4%) | 94 (0.4%) |
| burkina faso | 211 (2.1%) | 293 (2.3%) | 504 (2.2%) |
| cabo verde | 44 (0.4%) | 44 (0.3%) | 88 (0.4%) |
| cameroon | 319 (3.2%) | 640 (5.0%) | 959 (4.2%) |
| congo-brazzaville | 168 (1.7%) | 515 (4.0%) | 683 (3.0%) |
| côte d'ivoire | 445 (4.5%) | 318 (2.5%) | 763 (3.4%) |
| eswatini | 87 (0.9%) | 282 (2.2%) | 369 (1.6%) |
| ethiopia | 377 (3.8%) | 560 (4.4%) | 937 (4.1%) |
| gabon | 253 (2.6%) | 768 (6.0%) | 1021 (4.5%) |
| gambia | 176 (1.8%) | 407 (3.2%) | 583 (2.6%) |
| ghana | 305 (3.1%) | 648 (5.1%) | 953 (4.2%) |
| guinea | 158 (1.6%) | 102 (0.8%) | 260 (1.2%) |
| kenya | 1054 (10.7%) | 661 (5.2%) | 1715 (7.6%) |
| lesotho | 96 (1.0%) | 91 (0.7%) | 187 (0.8%) |
| liberia | 289 (2.9%) | 206 (1.6%) | 495 (2.2%) |
| madagascar | 370 (3.8%) | 537 (4.2%) | 907 (4.0%) |
| malawi | 366 (3.7%) | 458 (3.6%) | 824 (3.6%) |
| mali | 343 (3.5%) | 410 (3.2%) | 753 (3.3%) |
| mauritania | 140 (1.4%) | 226 (1.8%) | 366 (1.6%) |
| mauritius | 16 (0.2%) | 41 (0.3%) | 57 (0.3%) |
| morocco | 40 (0.4%) | 153 (1.2%) | 193 (0.9%) |
| mozambique | 161 (1.6%) | 151 (1.2%) | 312 (1.4%) |
| namibia | 323 (3.3%) | 519 (4.1%) | 842 (3.7%) |
| niger | 350 (3.6%) | 365 (2.9%) | 715 (3.2%) |
| nigeria | 490 (5.0%) | 541 (4.2%) | 1031 (4.6%) |
| são tomé and príncipe | 42 (0.4%) | 41 (0.3%) | 83 (0.4%) |
| senegal | 194 (2.0%) | 584 (4.6%) | 778 (3.4%) |
| seychelles | 33 (0.3%) | 148 (1.2%) | 181 (0.8%) |
| sierra leone | 259 (2.6%) | 184 (1.4%) | 443 (2.0%) |
| south africa | 128 (1.3%) | 486 (3.8%) | 614 (2.7%) |
| sudan | 323 (3.3%) | 356 (2.8%) | 679 (3.0%) |
| tanzania | 877 (8.9%) | 742 (5.8%) | 1619 (7.2%) |
| togo | 119 (1.2%) | 215 (1.7%) | 334 (1.5%) |
| tunisia | 29 (0.3%) | 142 (1.1%) | 171 (0.8%) |
| uganda | 542 (5.5%) | 130 (1.0%) | 672 (3.0%) |
| zambia | 115 (1.2%) | 255 (2.0%) | 370 (1.6%) |
| zimbabwe | 155 (1.6%) | 163 (1.3%) | 318 (1.4%) |
table1::table1(~ctry, afro)| Overall (N=53444) |
|
|---|---|
| ctry | |
| Angola | 1200 (2.2%) |
| Benin | 1200 (2.2%) |
| Botswana | 1200 (2.2%) |
| Burkina Faso | 1200 (2.2%) |
| Cabo Verde | 1199 (2.2%) |
| Cameroon | 1200 (2.2%) |
| Congo-Brazzaville | 1200 (2.2%) |
| Côte d'Ivoire | 1200 (2.2%) |
| Eswatini | 1200 (2.2%) |
| Ethiopia | 2400 (4.5%) |
| Gabon | 1200 (2.2%) |
| Gambia | 1200 (2.2%) |
| Ghana | 2369 (4.4%) |
| Guinea | 1200 (2.2%) |
| Kenya | 2400 (4.5%) |
| Lesotho | 1200 (2.2%) |
| Liberia | 1200 (2.2%) |
| Madagascar | 1200 (2.2%) |
| Malawi | 1200 (2.2%) |
| Mali | 1200 (2.2%) |
| Mauritania | 1200 (2.2%) |
| Mauritius | 1200 (2.2%) |
| Morocco | 1200 (2.2%) |
| Mozambique | 1120 (2.1%) |
| Namibia | 1200 (2.2%) |
| Niger | 1200 (2.2%) |
| Nigeria | 1600 (3.0%) |
| São Tomé and Príncipe | 1200 (2.2%) |
| Senegal | 1200 (2.2%) |
| Seychelles | 1176 (2.2%) |
| Sierra Leone | 1200 (2.2%) |
| South Africa | 1580 (3.0%) |
| Sudan | 1200 (2.2%) |
| Tanzania | 2400 (4.5%) |
| Togo | 1200 (2.2%) |
| Tunisia | 1200 (2.2%) |
| Uganda | 2400 (4.5%) |
| Zambia | 1200 (2.2%) |
| Zimbabwe | 1200 (2.2%) |
Firstly, we group language policy into three categories rather than four to allow the model to converge. Furthermore, we create a set of categorical variables for the respondents’ age, which will be used to compare categorical age and continuous age with regard to the interaction between age and language profile.
df_mlvl$lang_policy3 <- fct_collapse(
df_mlvl$indigenous_lang_integration,
"no integration" = "no integration",
"limited integration" = "limited integration",
"moderate or extensive" = c("moderate integration", "extensive integration")
)
# Sanity check
#table(df_mlvl$lang_policy3, df_mlvl$ctry) |> apply(1, \(x) sum(x > 0))
df_mlvl$age_group <- cut(df_mlvl$age,
breaks = c(17, 35, Inf),
labels = c("Youth", "Adult"))
df_mlvl <- df_mlvl %>%
mutate(
across(
where(is.numeric) & !all_of("service_utilization"),
~ as.numeric(scale(.))
)
)fit <- glm(hesitancy ~ 1, data = df_mlvl, family = binomial(link = "logit"))
null_model <- glmer(hesitancy ~ 1 + (1|ctry), data = df_mlvl,family = binomial(link = "logit"),
control = glmerControl(optimizer = "bobyqa"))
2*(logLik(null_model) - logLik(fit))'log Lik.' 1682.737 (df=2)icc(null_model)# Intraclass Correlation Coefficient
Adjusted ICC: 0.108
Unadjusted ICC: 0.108The intraclass correlation coefficient was 0.108, indicating that 10.8 percent of the total variance in hesitancy is attributable to differences between countries. This suggests that it makes sense for this study to apply hierarchical regression with country as a random intercept. Furthermore, the random effect suggests that there is variance in vaccine hesitancy within and between countries.
A likelihood ratio test comparing the null model to the full model showed a significant improvement in fit, χ²(2) = 1682.74, supporting the inclusion of a country level analysis.
model1 <- glmer(
hesitancy ~
language_profile +
# --- The Multilevel Component ---
(1 | ctry),
data = df_mlvl,
# Specify this is a logistic regression
family = binomial(link = "logit"),
control = glmerControl(optimizer = "bobyqa")
)anova(null_model, model1, test = "Chisq")Data: df_mlvl
Models:
null_model: hesitancy ~ 1 + (1 | ctry)
model1: hesitancy ~ language_profile + (1 | ctry)
npar AIC BIC logLik -2*log(L) Chisq Df Pr(>Chisq)
null_model 2 29260 29276 -14628 29256
model1 4 29185 29217 -14588 29177 79.344 2 < 2.2e-16 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1model2 <- glmer(
hesitancy ~
language_profile + avg_diff_model + avg_resp_model + service_utilization +
# --- The Multilevel Component ---
(1 | ctry),
data = df_mlvl,
# Specify this is a logistic regression
family = binomial(link = "logit"),
control = glmerControl(optimizer = "bobyqa")
)anova(model1, model2, test = "Chisq")Data: df_mlvl
Models:
model1: hesitancy ~ language_profile + (1 | ctry)
model2: hesitancy ~ language_profile + avg_diff_model + avg_resp_model + service_utilization + (1 | ctry)
npar AIC BIC logLik -2*log(L) Chisq Df Pr(>Chisq)
model1 4 29185 29217 -14588 29177
model2 7 29120 29176 -14553 29106 70.945 3 2.679e-15 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1model3 <- glmer(
hesitancy ~
language_profile + avg_diff_model + avg_resp_model + service_utilization +
Interpersonal_Index_dm + Institutional_Index_dm +
# --- The Multilevel Component ---
(1 | ctry),
data = df_mlvl,
# Specify this is a logistic regression
family = binomial(link = "logit"),
control = glmerControl(optimizer = "bobyqa")
)anova(model2, model3, test = "Chisq")Data: df_mlvl
Models:
model2: hesitancy ~ language_profile + avg_diff_model + avg_resp_model + service_utilization + (1 | ctry)
model3: hesitancy ~ language_profile + avg_diff_model + avg_resp_model + service_utilization + Interpersonal_Index_dm + Institutional_Index_dm + (1 | ctry)
npar AIC BIC logLik -2*log(L) Chisq Df Pr(>Chisq)
model2 7 29120 29176 -14553 29106
model3 9 28665 28737 -14323 28647 459.37 2 < 2.2e-16 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1model4 <- glmer(
hesitancy ~
language_profile + avg_diff_model + avg_resp_model + service_utilization +
Interpersonal_Index_dm + Institutional_Index_dm + civic_participation_score +
# --- The Multilevel Component ---
(1 | ctry),
data = df_mlvl,
# Specify this is a logistic regression
family = binomial(link = "logit"),
control = glmerControl(optimizer = "bobyqa")
)
anova(model3, model4, test = "Chisq")Data: df_mlvl
Models:
model3: hesitancy ~ language_profile + avg_diff_model + avg_resp_model + service_utilization + Interpersonal_Index_dm + Institutional_Index_dm + (1 | ctry)
model4: hesitancy ~ language_profile + avg_diff_model + avg_resp_model + service_utilization + Interpersonal_Index_dm + Institutional_Index_dm + civic_participation_score + (1 | ctry)
npar AIC BIC logLik -2*log(L) Chisq Df Pr(>Chisq)
model3 9 28665 28737 -14323 28647
model4 10 28606 28686 -14293 28586 61.022 1 5.645e-15 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1model5 <- glmer(
hesitancy ~
language_profile + avg_diff_model + avg_resp_model + service_utilization +
Interpersonal_Index_dm + Institutional_Index_dm + civic_participation_score + gndr + age + urban_rural + lvl_edu +
religion + lpi_dm +
# --- The Multilevel Component ---
(1 | ctry),
data = df_mlvl,
# Specify this is a logistic regression
family = binomial(link = "logit"),
control = glmerControl(optimizer = "bobyqa")
)
anova(model4, model5, test = "Chisq")Data: df_mlvl
Models:
model4: hesitancy ~ language_profile + avg_diff_model + avg_resp_model + service_utilization + Interpersonal_Index_dm + Institutional_Index_dm + civic_participation_score + (1 | ctry)
model5: hesitancy ~ language_profile + avg_diff_model + avg_resp_model + service_utilization + Interpersonal_Index_dm + Institutional_Index_dm + civic_participation_score + gndr + age + urban_rural + lvl_edu + religion + lpi_dm + (1 | ctry)
npar AIC BIC logLik -2*log(L) Chisq Df Pr(>Chisq)
model4 10 28606 28686 -14293 28586
model5 22 28537 28713 -14246 28493 93.117 12 1.227e-14 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1model_6 <- glmer(
hesitancy ~
language_profile + avg_diff_model + avg_resp_model + service_utilization +
Interpersonal_Index_dm + Institutional_Index_dm + civic_participation_score + gndr + age + urban_rural + lvl_edu +
religion + lpi_dm + lang_policy3 + ctry_Institutional_Index + ctry_Interpersonal_Index + ctry_lpi +
# --- The Multilevel Component ---
(1 | ctry),
data = df_mlvl,
# Specify this is a logistic regression
family = binomial(link = "logit"),
control = glmerControl(optimizer = "bobyqa")
)
saveRDS(model_6, "/Users/PaulMUTAMBA/Documents/Old Pc/KUL QASS/Master/R_data/model_6.rds")model_6 <- readRDS("/Users/PaulMUTAMBA/Documents/Old Pc/KUL QASS/Master/R_data/model_6.rds")
anova(model5, model_6)Data: df_mlvl
Models:
model5: hesitancy ~ language_profile + avg_diff_model + avg_resp_model + service_utilization + Interpersonal_Index_dm + Institutional_Index_dm + civic_participation_score + gndr + age + urban_rural + lvl_edu + religion + lpi_dm + (1 | ctry)
model_6: hesitancy ~ language_profile + avg_diff_model + avg_resp_model + service_utilization + Interpersonal_Index_dm + Institutional_Index_dm + civic_participation_score + gndr + age + urban_rural + lvl_edu + religion + lpi_dm + lang_policy3 + ctry_Institutional_Index + ctry_Interpersonal_Index + ctry_lpi + (1 | ctry)
npar AIC BIC logLik -2*log(L) Chisq Df Pr(>Chisq)
model5 22 28537 28713 -14246 28493
model_6 27 28539 28756 -14243 28485 7.3093 5 0.1986tab_model(model1, model2, model3,model4, model5,model_6)| hesitancy | hesitancy | hesitancy | hesitancy | hesitancy | hesitancy | |||||||||||||
|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|
| Predictors | Odds Ratios | CI | p | Odds Ratios | CI | p | Odds Ratios | CI | p | Odds Ratios | CI | p | Odds Ratios | CI | p | Odds Ratios | CI | p |
| (Intercept) | 1.95 | 1.55 – 2.44 | <0.001 | 1.78 | 1.38 – 2.30 | <0.001 | 1.65 | 1.28 – 2.14 | <0.001 | 1.58 | 1.23 – 2.04 | <0.001 | 1.67 | 1.25 – 2.24 | 0.001 | 1.83 | 1.31 – 2.57 | <0.001 |
| language profile [Consistent Indigenous] |
0.64 | 0.57 – 0.72 | <0.001 | 0.65 | 0.57 – 0.73 | <0.001 | 0.71 | 0.63 – 0.81 | <0.001 | 0.73 | 0.64 – 0.82 | <0.001 | 0.79 | 0.69 – 0.90 | <0.001 | 0.78 | 0.69 – 0.90 | <0.001 |
| language profile [Mixed/Code-Switcher] |
0.90 | 0.81 – 1.01 | 0.074 | 0.90 | 0.80 – 1.01 | 0.069 | 0.92 | 0.82 – 1.03 | 0.140 | 0.93 | 0.83 – 1.05 | 0.235 | 0.99 | 0.88 – 1.11 | 0.834 | 0.98 | 0.88 – 1.10 | 0.784 |
| avg diff model | 1.07 | 1.02 – 1.13 | 0.011 | 1.01 | 0.95 – 1.06 | 0.834 | 1.01 | 0.96 – 1.06 | 0.748 | 1.00 | 0.95 – 1.06 | 0.954 | 1.00 | 0.95 – 1.06 | 0.942 | |||
| avg resp model | 0.84 | 0.80 – 0.89 | <0.001 | 0.89 | 0.85 – 0.94 | <0.001 | 0.89 | 0.85 – 0.94 | <0.001 | 0.90 | 0.85 – 0.94 | <0.001 | 0.90 | 0.85 – 0.94 | <0.001 | |||
| service utilization | 1.12 | 0.95 – 1.31 | 0.173 | 1.14 | 0.97 – 1.34 | 0.116 | 1.17 | 1.00 – 1.38 | 0.051 | 1.18 | 1.01 – 1.39 | 0.043 | 1.08 | 1.00 – 1.16 | 0.044 | |||
| Interpersonal Index dm | 0.89 | 0.87 – 0.92 | <0.001 | 0.90 | 0.87 – 0.92 | <0.001 | 0.90 | 0.87 – 0.93 | <0.001 | 0.90 | 0.87 – 0.93 | <0.001 | ||||||
| Institutional Index dm | 0.77 | 0.75 – 0.79 | <0.001 | 0.77 | 0.75 – 0.79 | <0.001 | 0.78 | 0.75 – 0.80 | <0.001 | 0.78 | 0.75 – 0.80 | <0.001 | ||||||
| civic participation score | 0.89 | 0.86 – 0.92 | <0.001 | 0.89 | 0.86 – 0.92 | <0.001 | 0.89 | 0.86 – 0.92 | <0.001 | |||||||||
| gndr [woman] | 0.98 | 0.92 – 1.04 | 0.437 | 0.98 | 0.92 – 1.04 | 0.440 | ||||||||||||
| age | 1.03 | 1.00 – 1.06 | 0.038 | 1.03 | 1.00 – 1.06 | 0.042 | ||||||||||||
| urban rural [urban] | 1.23 | 1.15 – 1.31 | <0.001 | 1.23 | 1.15 – 1.31 | <0.001 | ||||||||||||
| lvl edu [Primary] | 0.98 | 0.90 – 1.07 | 0.674 | 0.98 | 0.89 – 1.07 | 0.638 | ||||||||||||
| lvl edu [Secondary] | 0.87 | 0.79 – 0.95 | 0.003 | 0.86 | 0.78 – 0.95 | 0.003 | ||||||||||||
| lvl edu [Higher Education] |
1.01 | 0.90 – 1.13 | 0.891 | 1.00 | 0.89 – 1.13 | 0.953 | ||||||||||||
| lvl edu [No Answer] | 0.91 | 0.48 – 1.74 | 0.781 | 0.91 | 0.47 – 1.73 | 0.764 | ||||||||||||
| religion [Muslim] | 0.73 | 0.63 – 0.85 | <0.001 | 0.73 | 0.63 – 0.85 | <0.001 | ||||||||||||
| religion [Protestant/Christian] |
0.92 | 0.81 – 1.05 | 0.216 | 0.92 | 0.81 – 1.05 | 0.211 | ||||||||||||
| religion [Other/No Answer] |
0.94 | 0.71 – 1.25 | 0.673 | 0.96 | 0.72 – 1.27 | 0.774 | ||||||||||||
| religion [Catholic] | 0.85 | 0.74 – 0.98 | 0.023 | 0.85 | 0.74 – 0.98 | 0.022 | ||||||||||||
| lpi dm | 1.02 | 0.99 – 1.05 | 0.246 | 1.02 | 0.99 – 1.05 | 0.248 | ||||||||||||
| lang policy3 [limited integration] |
0.92 | 0.58 – 1.45 | 0.710 | |||||||||||||||
| lang policy3 [moderate or extensive] |
1.16 | 0.72 – 1.86 | 0.541 | |||||||||||||||
| ctry Institutional Index | 0.83 | 0.62 – 1.12 | 0.220 | |||||||||||||||
| ctry Interpersonal Index | 1.08 | 0.87 – 1.34 | 0.500 | |||||||||||||||
| ctry lpi | 0.81 | 0.68 – 0.95 | 0.010 | |||||||||||||||
| Random Effects | ||||||||||||||||||
| σ2 | 3.29 | 3.29 | 3.29 | 3.29 | 3.29 | 3.29 | ||||||||||||
| τ00 | 0.40 ctry | 0.40 ctry | 0.42 ctry | 0.40 ctry | 0.40 ctry | 0.33 ctry | ||||||||||||
| ICC | 0.11 | 0.11 | 0.11 | 0.11 | 0.11 | 0.09 | ||||||||||||
| N | 39 ctry | 39 ctry | 39 ctry | 39 ctry | 39 ctry | 39 ctry | ||||||||||||
| Observations | 22598 | 22598 | 22598 | 22598 | 22598 | 22598 | ||||||||||||
| Marginal R2 / Conditional R2 | 0.009 / 0.117 | 0.014 / 0.122 | 0.035 / 0.143 | 0.039 / 0.144 | 0.046 / 0.150 | 0.059 / 0.146 | ||||||||||||
model_7 <- glmer(
hesitancy ~
language_profile + avg_diff_model + avg_resp_model + service_utilization +
Interpersonal_Index_dm + Institutional_Index_dm + civic_participation_score + gndr + age + urban_rural + lvl_edu +
religion + lpi_dm + lang_policy3 + ctry_Institutional_Index + ctry_Interpersonal_Index + ctry_lpi + language_profile:lang_policy3 +
# --- The Multilevel Component ---
(1 | ctry),
data = df_mlvl,
# Specify this is a logistic regression
family = binomial(link = "logit"),
control = glmerControl(optimizer = "bobyqa")
)
saveRDS(model_7, "/Users/PaulMUTAMBA/Documents/Old Pc/KUL QASS/Master/R_data/model_7.rds")model_7 <- readRDS("/Users/PaulMUTAMBA/Documents/Old Pc/KUL QASS/Master/R_data/model_7.rds")
anova(model_6, model_7)Data: df_mlvl
Models:
model_6: hesitancy ~ language_profile + avg_diff_model + avg_resp_model + service_utilization + Interpersonal_Index_dm + Institutional_Index_dm + civic_participation_score + gndr + age + urban_rural + lvl_edu + religion + lpi_dm + lang_policy3 + ctry_Institutional_Index + ctry_Interpersonal_Index + ctry_lpi + (1 | ctry)
model_7: hesitancy ~ language_profile + avg_diff_model + avg_resp_model + service_utilization + Interpersonal_Index_dm + Institutional_Index_dm + civic_participation_score + gndr + age + urban_rural + lvl_edu + religion + lpi_dm + lang_policy3 + ctry_Institutional_Index + ctry_Interpersonal_Index + ctry_lpi + language_profile:lang_policy3 + (1 | ctry)
npar AIC BIC logLik -2*log(L) Chisq Df Pr(>Chisq)
model_6 27 28539 28756 -14243 28485
model_7 31 28524 28773 -14231 28462 22.846 4 0.0001359 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1model_8 <- glmer(
hesitancy ~
language_profile + avg_diff_model + avg_resp_model + service_utilization +
Interpersonal_Index_dm + Institutional_Index_dm + civic_participation_score + gndr + age + urban_rural + lvl_edu +
religion + lpi_dm + lang_policy3 + ctry_Institutional_Index + ctry_Interpersonal_Index + ctry_lpi + language_profile:lang_policy3 + avg_resp_model:lang_policy3 +
# --- The Multilevel Component ---
(1 | ctry),
data = df_mlvl,
# Specify this is a logistic regression
family = binomial(link = "logit"),
control = glmerControl(optimizer = "bobyqa",
optCtrl = list(maxfun = 200000))
)
saveRDS(model_8, "/Users/PaulMUTAMBA/Documents/Old Pc/KUL QASS/Master/R_data/model_8.rds")model_8 <- readRDS("/Users/PaulMUTAMBA/Documents/Old Pc/KUL QASS/Master/R_data/model_8.rds")
anova(model_7, model_8)Data: df_mlvl
Models:
model_7: hesitancy ~ language_profile + avg_diff_model + avg_resp_model + service_utilization + Interpersonal_Index_dm + Institutional_Index_dm + civic_participation_score + gndr + age + urban_rural + lvl_edu + religion + lpi_dm + lang_policy3 + ctry_Institutional_Index + ctry_Interpersonal_Index + ctry_lpi + language_profile:lang_policy3 + (1 | ctry)
model_8: hesitancy ~ language_profile + avg_diff_model + avg_resp_model + service_utilization + Interpersonal_Index_dm + Institutional_Index_dm + civic_participation_score + gndr + age + urban_rural + lvl_edu + religion + lpi_dm + lang_policy3 + ctry_Institutional_Index + ctry_Interpersonal_Index + ctry_lpi + language_profile:lang_policy3 + avg_resp_model:lang_policy3 + (1 | ctry)
npar AIC BIC logLik -2*log(L) Chisq Df Pr(>Chisq)
model_7 31 28524 28773 -14231 28462
model_8 33 28527 28792 -14230 28461 1.4843 2 0.4761tab_model(model_6, model_7, model_8)| hesitancy | hesitancy | hesitancy | |||||||
|---|---|---|---|---|---|---|---|---|---|
| Predictors | Odds Ratios | CI | p | Odds Ratios | CI | p | Odds Ratios | CI | p |
| (Intercept) | 1.83 | 1.31 – 2.57 | <0.001 | 1.51 | 1.06 – 2.14 | 0.021 | 1.51 | 1.06 – 2.14 | 0.022 |
| language profile [Consistent Indigenous] |
0.78 | 0.69 – 0.90 | <0.001 | 1.05 | 0.87 – 1.27 | 0.632 | 1.05 | 0.87 – 1.27 | 0.627 |
| language profile [Mixed/Code-Switcher] |
0.98 | 0.88 – 1.10 | 0.784 | 0.99 | 0.86 – 1.15 | 0.922 | 0.99 | 0.86 – 1.15 | 0.931 |
| avg diff model | 1.00 | 0.95 – 1.06 | 0.942 | 1.00 | 0.95 – 1.06 | 0.916 | 1.00 | 0.95 – 1.06 | 0.969 |
| avg resp model | 0.90 | 0.85 – 0.94 | <0.001 | 0.90 | 0.85 – 0.95 | <0.001 | 0.92 | 0.86 – 0.98 | 0.011 |
| service utilization | 1.08 | 1.00 – 1.16 | 0.044 | 1.08 | 1.00 – 1.16 | 0.049 | 1.08 | 1.00 – 1.16 | 0.050 |
| Interpersonal Index dm | 0.90 | 0.87 – 0.93 | <0.001 | 0.90 | 0.88 – 0.93 | <0.001 | 0.90 | 0.88 – 0.93 | <0.001 |
| Institutional Index dm | 0.78 | 0.75 – 0.80 | <0.001 | 0.78 | 0.75 – 0.80 | <0.001 | 0.78 | 0.75 – 0.80 | <0.001 |
| civic participation score | 0.89 | 0.86 – 0.92 | <0.001 | 0.89 | 0.86 – 0.92 | <0.001 | 0.89 | 0.86 – 0.92 | <0.001 |
| gndr [woman] | 0.98 | 0.92 – 1.04 | 0.440 | 0.98 | 0.92 – 1.03 | 0.395 | 0.98 | 0.92 – 1.03 | 0.400 |
| age | 1.03 | 1.00 – 1.06 | 0.042 | 1.03 | 1.00 – 1.06 | 0.040 | 1.03 | 1.00 – 1.06 | 0.040 |
| urban rural [urban] | 1.23 | 1.15 – 1.31 | <0.001 | 1.23 | 1.15 – 1.31 | <0.001 | 1.23 | 1.15 – 1.31 | <0.001 |
| lvl edu [Primary] | 0.98 | 0.89 – 1.07 | 0.638 | 0.99 | 0.90 – 1.09 | 0.844 | 0.99 | 0.91 – 1.09 | 0.855 |
| lvl edu [Secondary] | 0.86 | 0.78 – 0.95 | 0.003 | 0.87 | 0.79 – 0.96 | 0.005 | 0.87 | 0.79 – 0.96 | 0.006 |
| lvl edu [Higher Education] |
1.00 | 0.89 – 1.13 | 0.953 | 1.01 | 0.90 – 1.13 | 0.900 | 1.01 | 0.90 – 1.13 | 0.899 |
| lvl edu [No Answer] | 0.91 | 0.47 – 1.73 | 0.764 | 0.94 | 0.49 – 1.80 | 0.843 | 0.94 | 0.49 – 1.79 | 0.841 |
| religion [Muslim] | 0.73 | 0.63 – 0.85 | <0.001 | 0.74 | 0.64 – 0.86 | <0.001 | 0.74 | 0.64 – 0.86 | <0.001 |
| religion [Protestant/Christian] |
0.92 | 0.81 – 1.05 | 0.211 | 0.93 | 0.81 – 1.05 | 0.237 | 0.93 | 0.81 – 1.05 | 0.240 |
| religion [Other/No Answer] |
0.96 | 0.72 – 1.27 | 0.774 | 0.95 | 0.72 – 1.26 | 0.736 | 0.95 | 0.72 – 1.26 | 0.730 |
| religion [Catholic] | 0.85 | 0.74 – 0.98 | 0.022 | 0.86 | 0.74 – 0.99 | 0.031 | 0.86 | 0.74 – 0.99 | 0.031 |
| lpi dm | 1.02 | 0.99 – 1.05 | 0.248 | 1.02 | 0.99 – 1.05 | 0.224 | 1.02 | 0.99 – 1.05 | 0.211 |
| lang policy3 [limited integration] |
0.92 | 0.58 – 1.45 | 0.710 | 1.38 | 0.83 – 2.31 | 0.219 | 1.38 | 0.82 – 2.31 | 0.220 |
| lang policy3 [moderate or extensive] |
1.16 | 0.72 – 1.86 | 0.541 | 1.36 | 0.76 – 2.43 | 0.303 | 1.36 | 0.76 – 2.44 | 0.300 |
| ctry Institutional Index | 0.83 | 0.62 – 1.12 | 0.220 | 0.82 | 0.61 – 1.09 | 0.178 | 0.82 | 0.61 – 1.10 | 0.179 |
| ctry Interpersonal Index | 1.08 | 0.87 – 1.34 | 0.500 | 1.08 | 0.87 – 1.35 | 0.474 | 1.08 | 0.87 – 1.35 | 0.476 |
| ctry lpi | 0.81 | 0.68 – 0.95 | 0.010 | 0.82 | 0.70 – 0.97 | 0.020 | 0.82 | 0.70 – 0.97 | 0.020 |
| language profile [Consistent Indigenous] × lang policy3 [limited integration] |
0.57 | 0.43 – 0.76 | <0.001 | 0.57 | 0.43 – 0.76 | <0.001 | |||
| language profile [Mixed/Code-Switcher] × lang policy3 [limited integration] |
0.88 | 0.67 – 1.16 | 0.373 | 0.88 | 0.67 – 1.16 | 0.369 | |||
| language profile [Consistent Indigenous] × lang policy3 [moderate or extensive] |
0.78 | 0.52 – 1.16 | 0.216 | 0.78 | 0.52 – 1.16 | 0.214 | |||
| language profile [Mixed/Code-Switcher] × lang policy3 [moderate or extensive] |
0.95 | 0.65 – 1.39 | 0.802 | 0.95 | 0.65 – 1.38 | 0.787 | |||
| avg resp model × lang policy3 [limited integration] |
0.98 | 0.91 – 1.05 | 0.509 | ||||||
| avg resp model × lang policy3 [moderate or extensive] |
0.96 | 0.89 – 1.03 | 0.224 | ||||||
| Random Effects | |||||||||
| σ2 | 3.29 | 3.29 | 3.29 | ||||||
| τ00 | 0.33 ctry | 0.33 ctry | 0.33 ctry | ||||||
| ICC | 0.09 | 0.09 | 0.09 | ||||||
| N | 39 ctry | 39 ctry | 39 ctry | ||||||
| Observations | 22598 | 22598 | 22598 | ||||||
| Marginal R2 / Conditional R2 | 0.059 / 0.146 | 0.060 / 0.145 | 0.060 / 0.145 | ||||||
tab_model(model_8, show.icc = TRUE, show.r2 = TRUE, show.aic = TRUE)| hesitancy | |||
|---|---|---|---|
| Predictors | Odds Ratios | CI | p |
| (Intercept) | 1.51 | 1.06 – 2.14 | 0.022 |
| language profile [Consistent Indigenous] |
1.05 | 0.87 – 1.27 | 0.627 |
| language profile [Mixed/Code-Switcher] |
0.99 | 0.86 – 1.15 | 0.931 |
| avg diff model | 1.00 | 0.95 – 1.06 | 0.969 |
| avg resp model | 0.92 | 0.86 – 0.98 | 0.011 |
| service utilization | 1.08 | 1.00 – 1.16 | 0.050 |
| Interpersonal Index dm | 0.90 | 0.88 – 0.93 | <0.001 |
| Institutional Index dm | 0.78 | 0.75 – 0.80 | <0.001 |
| civic participation score | 0.89 | 0.86 – 0.92 | <0.001 |
| gndr [woman] | 0.98 | 0.92 – 1.03 | 0.400 |
| age | 1.03 | 1.00 – 1.06 | 0.040 |
| urban rural [urban] | 1.23 | 1.15 – 1.31 | <0.001 |
| lvl edu [Primary] | 0.99 | 0.91 – 1.09 | 0.855 |
| lvl edu [Secondary] | 0.87 | 0.79 – 0.96 | 0.006 |
| lvl edu [Higher Education] |
1.01 | 0.90 – 1.13 | 0.899 |
| lvl edu [No Answer] | 0.94 | 0.49 – 1.79 | 0.841 |
| religion [Muslim] | 0.74 | 0.64 – 0.86 | <0.001 |
| religion [Protestant/Christian] |
0.93 | 0.81 – 1.05 | 0.240 |
| religion [Other/No Answer] |
0.95 | 0.72 – 1.26 | 0.730 |
| religion [Catholic] | 0.86 | 0.74 – 0.99 | 0.031 |
| lpi dm | 1.02 | 0.99 – 1.05 | 0.211 |
| lang policy3 [limited integration] |
1.38 | 0.82 – 2.31 | 0.220 |
| lang policy3 [moderate or extensive] |
1.36 | 0.76 – 2.44 | 0.300 |
| ctry Institutional Index | 0.82 | 0.61 – 1.10 | 0.179 |
| ctry Interpersonal Index | 1.08 | 0.87 – 1.35 | 0.476 |
| ctry lpi | 0.82 | 0.70 – 0.97 | 0.020 |
| language profile [Consistent Indigenous] × lang policy3 [limited integration] |
0.57 | 0.43 – 0.76 | <0.001 |
| language profile [Mixed/Code-Switcher] × lang policy3 [limited integration] |
0.88 | 0.67 – 1.16 | 0.369 |
| language profile [Consistent Indigenous] × lang policy3 [moderate or extensive] |
0.78 | 0.52 – 1.16 | 0.214 |
| language profile [Mixed/Code-Switcher] × lang policy3 [moderate or extensive] |
0.95 | 0.65 – 1.38 | 0.787 |
| avg resp model × lang policy3 [limited integration] |
0.98 | 0.91 – 1.05 | 0.509 |
| avg resp model × lang policy3 [moderate or extensive] |
0.96 | 0.89 – 1.03 | 0.224 |
| Random Effects | |||
| σ2 | 3.29 | ||
| τ00 ctry | 0.33 | ||
| ICC | 0.09 | ||
| N ctry | 39 | ||
| Observations | 22598 | ||
| Marginal R2 / Conditional R2 | 0.060 / 0.145 | ||
| AIC | 28527.001 | ||
Below are the countries effect
country_effects <- lme4::ranef(model_8, condVar = TRUE)
##-----Plot the caterpillar plots----##
lattice::dotplot(country_effects)$ctry
pred_lang <- ggpredict(model_8, terms = c("language_profile", "lang_policy3"))
plot(pred_lang) +
labs(title = NULL,
x = "Language profile",
y = "Predicted probability of hesitancy",
colour = "National language policy")
model_9 <- glmer(
hesitancy ~
(language_profile + avg_diff_model + avg_resp_model + service_utilization +
Interpersonal_Index_dm + Institutional_Index_dm) * age_group + civic_participation_score + gndr + urban_rural + lvl_edu +
religion + lpi_dm + lang_policy3 + ctry_Institutional_Index + ctry_Interpersonal_Index + ctry_lpi + language_profile:lang_policy3 + avg_resp_model:lang_policy3 +
# --- The Multilevel Component ---
(1 | ctry),
data = df_mlvl,
# Specify this is a logistic regression
family = binomial(link = "logit"),
control = glmerControl(optimizer = "bobyqa",
optCtrl = list(maxfun = 200000))
)
saveRDS(model_9, "/Users/PaulMUTAMBA/Documents/Old Pc/KUL QASS/Master/R_data/model_9.rds")model_9 <- readRDS("/Users/PaulMUTAMBA/Documents/Old Pc/KUL QASS/Master/R_data/model_9.rds")
pred_resp_age <- ggpredict(model_9, terms = c("avg_resp_model", "age_group"))
plot(pred_resp_age) +
labs(title = NULL,
x = "Perceived institutional respect (standardized)",
y = "Predicted probability of hesitancy",
colour = "Age group")
model_9 <- readRDS("/Users/PaulMUTAMBA/Documents/Old Pc/KUL QASS/Master/R_data/model_9.rds")
tab_model(model_9)| hesitancy | |||
|---|---|---|---|
| Predictors | Odds Ratios | CI | p |
| (Intercept) | 1.43 | 1.00 – 2.05 | 0.047 |
| language profile [Consistent Indigenous] |
1.09 | 0.89 – 1.34 | 0.398 |
| language profile [Mixed/Code-Switcher] |
1.00 | 0.85 – 1.18 | 0.981 |
| avg diff model | 1.00 | 0.94 – 1.08 | 0.912 |
| avg resp model | 0.87 | 0.81 – 0.94 | 0.001 |
| service utilization | 1.13 | 1.03 – 1.24 | 0.012 |
| Interpersonal Index dm | 0.91 | 0.87 – 0.94 | <0.001 |
| Institutional Index dm | 0.77 | 0.74 – 0.80 | <0.001 |
| age group [Adult] | 1.14 | 0.96 – 1.35 | 0.124 |
| civic participation score | 0.89 | 0.86 – 0.92 | <0.001 |
| gndr [woman] | 0.97 | 0.92 – 1.03 | 0.376 |
| urban rural [urban] | 1.23 | 1.15 – 1.31 | <0.001 |
| lvl edu [Primary] | 0.99 | 0.90 – 1.08 | 0.792 |
| lvl edu [Secondary] | 0.87 | 0.79 – 0.96 | 0.004 |
| lvl edu [Higher Education] |
1.00 | 0.89 – 1.13 | 0.967 |
| lvl edu [No Answer] | 0.94 | 0.49 – 1.80 | 0.847 |
| religion [Muslim] | 0.74 | 0.64 – 0.86 | <0.001 |
| religion [Protestant/Christian] |
0.93 | 0.81 – 1.05 | 0.243 |
| religion [Other/No Answer] |
0.95 | 0.71 – 1.26 | 0.706 |
| religion [Catholic] | 0.86 | 0.74 – 0.99 | 0.032 |
| lpi dm | 1.02 | 0.99 – 1.05 | 0.195 |
| lang policy3 [limited integration] |
1.38 | 0.83 – 2.32 | 0.218 |
| lang policy3 [moderate or extensive] |
1.37 | 0.76 – 2.46 | 0.293 |
| ctry Institutional Index | 0.82 | 0.61 – 1.10 | 0.183 |
| ctry Interpersonal Index | 1.08 | 0.87 – 1.35 | 0.477 |
| ctry lpi | 0.82 | 0.70 – 0.97 | 0.020 |
| language profile [Consistent Indigenous] × age group [Adult] |
0.90 | 0.75 – 1.08 | 0.267 |
| language profile [Mixed/Code-Switcher] × age group [Adult] |
0.99 | 0.80 – 1.22 | 0.911 |
| avg diff model × age group [Adult] |
1.00 | 0.89 – 1.11 | 0.931 |
| avg resp model × age group [Adult] |
1.13 | 1.02 – 1.25 | 0.016 |
| service utilization × age group [Adult] |
0.89 | 0.77 – 1.02 | 0.103 |
| Interpersonal Index dm × age group [Adult] |
0.99 | 0.93 – 1.05 | 0.760 |
| Institutional Index dm × age group [Adult] |
1.01 | 0.96 – 1.08 | 0.660 |
| language profile [Consistent Indigenous] × lang policy3 [limited integration] |
0.57 | 0.43 – 0.76 | <0.001 |
| language profile [Mixed/Code-Switcher] × lang policy3 [limited integration] |
0.88 | 0.67 – 1.16 | 0.374 |
| language profile [Consistent Indigenous] × lang policy3 [moderate or extensive] |
0.77 | 0.52 – 1.15 | 0.204 |
| language profile [Mixed/Code-Switcher] × lang policy3 [moderate or extensive] |
0.95 | 0.65 – 1.38 | 0.786 |
| avg resp model × lang policy3 [limited integration] |
0.98 | 0.91 – 1.05 | 0.510 |
| avg resp model × lang policy3 [moderate or extensive] |
0.96 | 0.89 – 1.03 | 0.225 |
| Random Effects | |||
| σ2 | 3.29 | ||
| τ00 ctry | 0.33 | ||
| ICC | 0.09 | ||
| N ctry | 39 | ||
| Observations | 22598 | ||
| Marginal R2 / Conditional R2 | 0.060 / 0.146 | ||
anova(model_8, model_9)Data: df_mlvl
Models:
model_8: hesitancy ~ language_profile + avg_diff_model + avg_resp_model + service_utilization + Interpersonal_Index_dm + Institutional_Index_dm + civic_participation_score + gndr + age + urban_rural + lvl_edu + religion + lpi_dm + lang_policy3 + ctry_Institutional_Index + ctry_Interpersonal_Index + ctry_lpi + language_profile:lang_policy3 + avg_resp_model:lang_policy3 + (1 | ctry)
model_9: hesitancy ~ (language_profile + avg_diff_model + avg_resp_model + service_utilization + Interpersonal_Index_dm + Institutional_Index_dm) * age_group + civic_participation_score + gndr + urban_rural + lvl_edu + religion + lpi_dm + lang_policy3 + ctry_Institutional_Index + ctry_Interpersonal_Index + ctry_lpi + language_profile:lang_policy3 + avg_resp_model:lang_policy3 + (1 | ctry)
npar AIC BIC logLik -2*log(L) Chisq Df Pr(>Chisq)
model_8 33 28527 28792 -14230 28461
model_9 40 28532 28853 -14226 28452 9.1399 7 0.2428tab_model(null_model, model5, model_6, model_8, model_9,
dv.labels = c("Null", "Individual + Controls",
"+ Country level", "+ Cross-level interactions",
"+ Age moderation"),
show.icc = TRUE, show.r2 = TRUE, show.aic = TRUE,
transform = "exp")| Null | Individual + Controls | + Country level | + Cross-level interactions | + Age moderation | |||||||||||
|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|
| Predictors | Odds Ratios | CI | p | Odds Ratios | CI | p | Odds Ratios | CI | p | Odds Ratios | CI | p | Odds Ratios | CI | p |
| (Intercept) | 1.41 | 1.15 – 1.73 | 0.001 | 1.67 | 1.25 – 2.24 | 0.001 | 1.83 | 1.31 – 2.57 | <0.001 | 1.51 | 1.06 – 2.14 | 0.022 | 1.43 | 1.00 – 2.05 | 0.047 |
| language profile [Consistent Indigenous] |
0.79 | 0.69 – 0.90 | <0.001 | 0.78 | 0.69 – 0.90 | <0.001 | 1.05 | 0.87 – 1.27 | 0.627 | 1.09 | 0.89 – 1.34 | 0.398 | |||
| language profile [Mixed/Code-Switcher] |
0.99 | 0.88 – 1.11 | 0.834 | 0.98 | 0.88 – 1.10 | 0.784 | 0.99 | 0.86 – 1.15 | 0.931 | 1.00 | 0.85 – 1.18 | 0.981 | |||
| avg diff model | 1.00 | 0.95 – 1.06 | 0.954 | 1.00 | 0.95 – 1.06 | 0.942 | 1.00 | 0.95 – 1.06 | 0.969 | 1.00 | 0.94 – 1.08 | 0.912 | |||
| avg resp model | 0.90 | 0.85 – 0.94 | <0.001 | 0.90 | 0.85 – 0.94 | <0.001 | 0.92 | 0.86 – 0.98 | 0.011 | 0.87 | 0.81 – 0.94 | 0.001 | |||
| service utilization | 1.18 | 1.01 – 1.39 | 0.043 | 1.08 | 1.00 – 1.16 | 0.044 | 1.08 | 1.00 – 1.16 | 0.050 | 1.13 | 1.03 – 1.24 | 0.012 | |||
| Interpersonal Index dm | 0.90 | 0.87 – 0.93 | <0.001 | 0.90 | 0.87 – 0.93 | <0.001 | 0.90 | 0.88 – 0.93 | <0.001 | 0.91 | 0.87 – 0.94 | <0.001 | |||
| Institutional Index dm | 0.78 | 0.75 – 0.80 | <0.001 | 0.78 | 0.75 – 0.80 | <0.001 | 0.78 | 0.75 – 0.80 | <0.001 | 0.77 | 0.74 – 0.80 | <0.001 | |||
| civic participation score | 0.89 | 0.86 – 0.92 | <0.001 | 0.89 | 0.86 – 0.92 | <0.001 | 0.89 | 0.86 – 0.92 | <0.001 | 0.89 | 0.86 – 0.92 | <0.001 | |||
| gndr [woman] | 0.98 | 0.92 – 1.04 | 0.437 | 0.98 | 0.92 – 1.04 | 0.440 | 0.98 | 0.92 – 1.03 | 0.400 | 0.97 | 0.92 – 1.03 | 0.376 | |||
| age | 1.03 | 1.00 – 1.06 | 0.038 | 1.03 | 1.00 – 1.06 | 0.042 | 1.03 | 1.00 – 1.06 | 0.040 | ||||||
| urban rural [urban] | 1.23 | 1.15 – 1.31 | <0.001 | 1.23 | 1.15 – 1.31 | <0.001 | 1.23 | 1.15 – 1.31 | <0.001 | 1.23 | 1.15 – 1.31 | <0.001 | |||
| lvl edu [Primary] | 0.98 | 0.90 – 1.07 | 0.674 | 0.98 | 0.89 – 1.07 | 0.638 | 0.99 | 0.91 – 1.09 | 0.855 | 0.99 | 0.90 – 1.08 | 0.792 | |||
| lvl edu [Secondary] | 0.87 | 0.79 – 0.95 | 0.003 | 0.86 | 0.78 – 0.95 | 0.003 | 0.87 | 0.79 – 0.96 | 0.006 | 0.87 | 0.79 – 0.96 | 0.004 | |||
| lvl edu [Higher Education] |
1.01 | 0.90 – 1.13 | 0.891 | 1.00 | 0.89 – 1.13 | 0.953 | 1.01 | 0.90 – 1.13 | 0.899 | 1.00 | 0.89 – 1.13 | 0.967 | |||
| lvl edu [No Answer] | 0.91 | 0.48 – 1.74 | 0.781 | 0.91 | 0.47 – 1.73 | 0.764 | 0.94 | 0.49 – 1.79 | 0.841 | 0.94 | 0.49 – 1.80 | 0.847 | |||
| religion [Muslim] | 0.73 | 0.63 – 0.85 | <0.001 | 0.73 | 0.63 – 0.85 | <0.001 | 0.74 | 0.64 – 0.86 | <0.001 | 0.74 | 0.64 – 0.86 | <0.001 | |||
| religion [Protestant/Christian] |
0.92 | 0.81 – 1.05 | 0.216 | 0.92 | 0.81 – 1.05 | 0.211 | 0.93 | 0.81 – 1.05 | 0.240 | 0.93 | 0.81 – 1.05 | 0.243 | |||
| religion [Other/No Answer] |
0.94 | 0.71 – 1.25 | 0.673 | 0.96 | 0.72 – 1.27 | 0.774 | 0.95 | 0.72 – 1.26 | 0.730 | 0.95 | 0.71 – 1.26 | 0.706 | |||
| religion [Catholic] | 0.85 | 0.74 – 0.98 | 0.023 | 0.85 | 0.74 – 0.98 | 0.022 | 0.86 | 0.74 – 0.99 | 0.031 | 0.86 | 0.74 – 0.99 | 0.032 | |||
| lpi dm | 1.02 | 0.99 – 1.05 | 0.246 | 1.02 | 0.99 – 1.05 | 0.248 | 1.02 | 0.99 – 1.05 | 0.211 | 1.02 | 0.99 – 1.05 | 0.195 | |||
| lang policy3 [limited integration] |
0.92 | 0.58 – 1.45 | 0.710 | 1.38 | 0.82 – 2.31 | 0.220 | 1.38 | 0.83 – 2.32 | 0.218 | ||||||
| lang policy3 [moderate or extensive] |
1.16 | 0.72 – 1.86 | 0.541 | 1.36 | 0.76 – 2.44 | 0.300 | 1.37 | 0.76 – 2.46 | 0.293 | ||||||
| ctry Institutional Index | 0.83 | 0.62 – 1.12 | 0.220 | 0.82 | 0.61 – 1.10 | 0.179 | 0.82 | 0.61 – 1.10 | 0.183 | ||||||
| ctry Interpersonal Index | 1.08 | 0.87 – 1.34 | 0.500 | 1.08 | 0.87 – 1.35 | 0.476 | 1.08 | 0.87 – 1.35 | 0.477 | ||||||
| ctry lpi | 0.81 | 0.68 – 0.95 | 0.010 | 0.82 | 0.70 – 0.97 | 0.020 | 0.82 | 0.70 – 0.97 | 0.020 | ||||||
| language profile [Consistent Indigenous] × lang policy3 [limited integration] |
0.57 | 0.43 – 0.76 | <0.001 | 0.57 | 0.43 – 0.76 | <0.001 | |||||||||
| language profile [Mixed/Code-Switcher] × lang policy3 [limited integration] |
0.88 | 0.67 – 1.16 | 0.369 | 0.88 | 0.67 – 1.16 | 0.374 | |||||||||
| language profile [Consistent Indigenous] × lang policy3 [moderate or extensive] |
0.78 | 0.52 – 1.16 | 0.214 | 0.77 | 0.52 – 1.15 | 0.204 | |||||||||
| language profile [Mixed/Code-Switcher] × lang policy3 [moderate or extensive] |
0.95 | 0.65 – 1.38 | 0.787 | 0.95 | 0.65 – 1.38 | 0.786 | |||||||||
| avg resp model × lang policy3 [limited integration] |
0.98 | 0.91 – 1.05 | 0.509 | 0.98 | 0.91 – 1.05 | 0.510 | |||||||||
| avg resp model × lang policy3 [moderate or extensive] |
0.96 | 0.89 – 1.03 | 0.224 | 0.96 | 0.89 – 1.03 | 0.225 | |||||||||
| age group [Adult] | 1.14 | 0.96 – 1.35 | 0.124 | ||||||||||||
| language profile [Consistent Indigenous] × age group [Adult] |
0.90 | 0.75 – 1.08 | 0.267 | ||||||||||||
| language profile [Mixed/Code-Switcher] × age group [Adult] |
0.99 | 0.80 – 1.22 | 0.911 | ||||||||||||
| avg diff model × age group [Adult] |
1.00 | 0.89 – 1.11 | 0.931 | ||||||||||||
| avg resp model × age group [Adult] |
1.13 | 1.02 – 1.25 | 0.016 | ||||||||||||
| service utilization × age group [Adult] |
0.89 | 0.77 – 1.02 | 0.103 | ||||||||||||
| Interpersonal Index dm × age group [Adult] |
0.99 | 0.93 – 1.05 | 0.760 | ||||||||||||
| Institutional Index dm × age group [Adult] |
1.01 | 0.96 – 1.08 | 0.660 | ||||||||||||
| Random Effects | |||||||||||||||
| σ2 | 3.29 | 3.29 | 3.29 | 3.29 | 3.29 | ||||||||||
| τ00 | 0.40 ctry | 0.40 ctry | 0.33 ctry | 0.33 ctry | 0.33 ctry | ||||||||||
| ICC | 0.11 | 0.11 | 0.09 | 0.09 | 0.09 | ||||||||||
| N | 39 ctry | 39 ctry | 39 ctry | 39 ctry | 39 ctry | ||||||||||
| Observations | 22598 | 22598 | 22598 | 22598 | 22598 | ||||||||||
| Marginal R2 / Conditional R2 | 0.000 / 0.108 | 0.046 / 0.150 | 0.059 / 0.146 | 0.060 / 0.145 | 0.060 / 0.146 | ||||||||||
| AIC | 29260.438 | 28536.640 | 28539.331 | 28527.001 | 28531.861 | ||||||||||
pred_age <- ggpredict(model_9, terms = c("language_profile", "age_group"))
plot(pred_age)Ignoring unknown labels:
• linetype : "age_group"
• shape : "age_group"
model_10 <- glmer(
hesitancy ~
(language_profile + avg_diff_model + avg_resp_model + service_utilization +
Interpersonal_Index_dm + Institutional_Index_dm) * age + civic_participation_score + gndr + urban_rural + lvl_edu +
religion + lpi_dm + lang_policy3 + ctry_Institutional_Index + ctry_Interpersonal_Index + ctry_lpi + language_profile:lang_policy3 + avg_resp_model:lang_policy3 +
# --- The Multilevel Component ---
(1 | ctry),
data = df_mlvl,
# Specify this is a logistic regression
family = binomial(link = "logit"),
control = glmerControl(optimizer = "bobyqa",
optCtrl = list(maxfun = 200000))
)
saveRDS(model_10, "/Users/PaulMUTAMBA/Documents/Old Pc/KUL QASS/Master/R_data/model_10.rds")model_10 <- readRDS("/Users/PaulMUTAMBA/Documents/Old Pc/KUL QASS/Master/R_data/model_10.rds")
tab_model(model_10)| hesitancy | |||
|---|---|---|---|
| Predictors | Odds Ratios | CI | p |
| (Intercept) | 1.37 | 0.95 – 1.98 | 0.095 |
| language profile [Consistent Indigenous] |
1.03 | 0.85 – 1.25 | 0.728 |
| language profile [Mixed/Code-Switcher] |
0.98 | 0.85 – 1.14 | 0.795 |
| avg diff model | 1.00 | 0.95 – 1.06 | 0.970 |
| avg resp model | 0.92 | 0.86 – 0.98 | 0.012 |
| service utilization | 1.17 | 1.00 – 1.38 | 0.053 |
| Interpersonal Index dm | 0.90 | 0.88 – 0.93 | <0.001 |
| Institutional Index dm | 0.78 | 0.75 – 0.80 | <0.001 |
| age | 1.21 | 1.05 – 1.41 | 0.010 |
| civic participation score | 0.89 | 0.86 – 0.92 | <0.001 |
| gndr [woman] | 0.97 | 0.92 – 1.03 | 0.389 |
| urban rural [urban] | 1.23 | 1.15 – 1.31 | <0.001 |
| lvl edu [Primary] | 0.99 | 0.90 – 1.08 | 0.812 |
| lvl edu [Secondary] | 0.87 | 0.79 – 0.96 | 0.005 |
| lvl edu [Higher Education] |
1.00 | 0.89 – 1.13 | 0.974 |
| lvl edu [No Answer] | 0.93 | 0.49 – 1.78 | 0.827 |
| religion [Muslim] | 0.74 | 0.64 – 0.85 | <0.001 |
| religion [Protestant/Christian] |
0.92 | 0.81 – 1.05 | 0.230 |
| religion [Other/No Answer] |
0.94 | 0.71 – 1.25 | 0.693 |
| religion [Catholic] | 0.85 | 0.74 – 0.98 | 0.027 |
| lpi dm | 1.02 | 0.99 – 1.05 | 0.192 |
| lang policy3 [limited integration] |
1.38 | 0.83 – 2.32 | 0.218 |
| lang policy3 [moderate or extensive] |
1.39 | 0.78 – 2.50 | 0.268 |
| ctry Institutional Index | 0.82 | 0.61 – 1.10 | 0.183 |
| ctry Interpersonal Index | 1.08 | 0.87 – 1.35 | 0.478 |
| ctry lpi | 0.82 | 0.70 – 0.97 | 0.020 |
| language profile [Consistent Indigenous] × age |
0.90 | 0.82 – 1.00 | 0.039 |
| language profile [Mixed/Code-Switcher] × age |
0.93 | 0.83 – 1.04 | 0.205 |
| avg diff model × age | 0.99 | 0.94 – 1.04 | 0.612 |
| avg resp model × age | 1.05 | 1.00 – 1.11 | 0.045 |
| service utilization × age | 0.90 | 0.76 – 1.05 | 0.182 |
| Interpersonal Index dm × age |
0.99 | 0.96 – 1.02 | 0.346 |
| Institutional Index dm × age |
1.01 | 0.98 – 1.04 | 0.572 |
| language profile [Consistent Indigenous] × lang policy3 [limited integration] |
0.57 | 0.43 – 0.77 | <0.001 |
| language profile [Mixed/Code-Switcher] × lang policy3 [limited integration] |
0.88 | 0.67 – 1.16 | 0.372 |
| language profile [Consistent Indigenous] × lang policy3 [moderate or extensive] |
0.76 | 0.51 – 1.13 | 0.176 |
| language profile [Mixed/Code-Switcher] × lang policy3 [moderate or extensive] |
0.93 | 0.64 – 1.36 | 0.717 |
| avg resp model × lang policy3 [limited integration] |
0.98 | 0.91 – 1.05 | 0.508 |
| avg resp model × lang policy3 [moderate or extensive] |
0.96 | 0.89 – 1.03 | 0.224 |
| Random Effects | |||
| σ2 | 3.29 | ||
| τ00 ctry | 0.33 | ||
| ICC | 0.09 | ||
| N ctry | 39 | ||
| Observations | 22598 | ||
| Marginal R2 / Conditional R2 | 0.060 / 0.146 | ||
summary(model_10)Generalized linear mixed model fit by maximum likelihood (Laplace
Approximation) [glmerMod]
Family: binomial ( logit )
Formula: hesitancy ~ (language_profile + avg_diff_model + avg_resp_model +
service_utilization + Interpersonal_Index_dm + Institutional_Index_dm) *
age + civic_participation_score + gndr + urban_rural + lvl_edu +
religion + lpi_dm + lang_policy3 + ctry_Institutional_Index +
ctry_Interpersonal_Index + ctry_lpi + language_profile:lang_policy3 +
avg_resp_model:lang_policy3 + (1 | ctry)
Data: df_mlvl
Control: glmerControl(optimizer = "bobyqa", optCtrl = list(maxfun = 2e+05))
AIC BIC logLik -2*log(L) df.resid
28529.2 28850.2 -14224.6 28449.2 22558
Scaled residuals:
Min 1Q Median 3Q Max
-3.2153 -0.9356 0.5103 0.8181 2.9908
Random effects:
Groups Name Variance Std.Dev.
ctry (Intercept) 0.3299 0.5744
Number of obs: 22598, groups: ctry, 39
Fixed effects:
Estimate
(Intercept) 0.313382
language_profileConsistent Indigenous 0.034040
language_profileMixed/Code-Switcher -0.019527
avg_diff_model 0.001033
avg_resp_model -0.082795
service_utilization 0.160388
Interpersonal_Index_dm -0.101640
Institutional_Index_dm -0.253427
age 0.193820
civic_participation_score -0.118365
gndrwoman -0.025449
urban_ruralurban 0.205953
lvl_eduPrimary -0.011078
lvl_eduSecondary -0.138693
lvl_eduHigher Education 0.001985
lvl_eduNo Answer -0.072660
religionMuslim -0.303395
religionProtestant/Christian -0.078810
religionOther/No Answer -0.056974
religionCatholic -0.159071
lpi_dm 0.019295
lang_policy3limited integration 0.324504
lang_policy3moderate or extensive 0.330996
ctry_Institutional_Index -0.196627
ctry_Interpersonal_Index 0.079230
ctry_lpi -0.193369
language_profileConsistent Indigenous:age -0.102566
language_profileMixed/Code-Switcher:age -0.073992
avg_diff_model:age -0.013622
avg_resp_model:age 0.050492
service_utilization:age -0.108846
Interpersonal_Index_dm:age -0.014273
Institutional_Index_dm:age 0.008333
language_profileConsistent Indigenous:lang_policy3limited integration -0.560950
language_profileMixed/Code-Switcher:lang_policy3limited integration -0.123835
language_profileConsistent Indigenous:lang_policy3moderate or extensive -0.274198
language_profileMixed/Code-Switcher:lang_policy3moderate or extensive -0.070010
avg_resp_model:lang_policy3limited integration -0.023561
avg_resp_model:lang_policy3moderate or extensive -0.044115
Std. Error
(Intercept) 0.187611
language_profileConsistent Indigenous 0.097815
language_profileMixed/Code-Switcher 0.075222
avg_diff_model 0.027444
avg_resp_model 0.033104
service_utilization 0.082747
Interpersonal_Index_dm 0.015006
Institutional_Index_dm 0.015151
age 0.074946
civic_participation_score 0.015516
gndrwoman 0.029566
urban_ruralurban 0.032706
lvl_eduPrimary 0.046648
lvl_eduSecondary 0.049121
lvl_eduHigher Education 0.060474
lvl_eduNo Answer 0.332085
religionMuslim 0.074172
religionProtestant/Christian 0.065610
religionOther/No Answer 0.144152
religionCatholic 0.071903
lpi_dm 0.014805
lang_policy3limited integration 0.263585
lang_policy3moderate or extensive 0.298610
ctry_Institutional_Index 0.147807
ctry_Interpersonal_Index 0.111553
ctry_lpi 0.083117
language_profileConsistent Indigenous:age 0.049773
language_profileMixed/Code-Switcher:age 0.058393
avg_diff_model:age 0.026893
avg_resp_model:age 0.025222
service_utilization:age 0.081582
Interpersonal_Index_dm:age 0.015149
Institutional_Index_dm:age 0.014760
language_profileConsistent Indigenous:lang_policy3limited integration 0.149800
language_profileMixed/Code-Switcher:lang_policy3limited integration 0.138849
language_profileConsistent Indigenous:lang_policy3moderate or extensive 0.202466
language_profileMixed/Code-Switcher:lang_policy3moderate or extensive 0.192940
avg_resp_model:lang_policy3limited integration 0.035575
avg_resp_model:lang_policy3moderate or extensive 0.036289
z value
(Intercept) 1.670
language_profileConsistent Indigenous 0.348
language_profileMixed/Code-Switcher -0.260
avg_diff_model 0.038
avg_resp_model -2.501
service_utilization 1.938
Interpersonal_Index_dm -6.773
Institutional_Index_dm -16.726
age 2.586
civic_participation_score -7.629
gndrwoman -0.861
urban_ruralurban 6.297
lvl_eduPrimary -0.237
lvl_eduSecondary -2.824
lvl_eduHigher Education 0.033
lvl_eduNo Answer -0.219
religionMuslim -4.090
religionProtestant/Christian -1.201
religionOther/No Answer -0.395
religionCatholic -2.212
lpi_dm 1.303
lang_policy3limited integration 1.231
lang_policy3moderate or extensive 1.108
ctry_Institutional_Index -1.330
ctry_Interpersonal_Index 0.710
ctry_lpi -2.326
language_profileConsistent Indigenous:age -2.061
language_profileMixed/Code-Switcher:age -1.267
avg_diff_model:age -0.507
avg_resp_model:age 2.002
service_utilization:age -1.334
Interpersonal_Index_dm:age -0.942
Institutional_Index_dm:age 0.565
language_profileConsistent Indigenous:lang_policy3limited integration -3.745
language_profileMixed/Code-Switcher:lang_policy3limited integration -0.892
language_profileConsistent Indigenous:lang_policy3moderate or extensive -1.354
language_profileMixed/Code-Switcher:lang_policy3moderate or extensive -0.363
avg_resp_model:lang_policy3limited integration -0.662
avg_resp_model:lang_policy3moderate or extensive -1.216
Pr(>|z|)
(Intercept) 0.094844
language_profileConsistent Indigenous 0.727834
language_profileMixed/Code-Switcher 0.795176
avg_diff_model 0.969973
avg_resp_model 0.012381
service_utilization 0.052587
Interpersonal_Index_dm 1.26e-11
Institutional_Index_dm < 2e-16
age 0.009706
civic_participation_score 2.37e-14
gndrwoman 0.389367
urban_ruralurban 3.03e-10
lvl_eduPrimary 0.812281
lvl_eduSecondary 0.004750
lvl_eduHigher Education 0.973816
lvl_eduNo Answer 0.826805
religionMuslim 4.31e-05
religionProtestant/Christian 0.229678
religionOther/No Answer 0.692669
religionCatholic 0.026946
lpi_dm 0.192460
lang_policy3limited integration 0.218278
lang_policy3moderate or extensive 0.267666
ctry_Institutional_Index 0.183422
ctry_Interpersonal_Index 0.477551
ctry_lpi 0.019994
language_profileConsistent Indigenous:age 0.039336
language_profileMixed/Code-Switcher:age 0.205100
avg_diff_model:age 0.612489
avg_resp_model:age 0.045295
service_utilization:age 0.182142
Interpersonal_Index_dm:age 0.346114
Institutional_Index_dm:age 0.572350
language_profileConsistent Indigenous:lang_policy3limited integration 0.000181
language_profileMixed/Code-Switcher:lang_policy3limited integration 0.372467
language_profileConsistent Indigenous:lang_policy3moderate or extensive 0.175644
language_profileMixed/Code-Switcher:lang_policy3moderate or extensive 0.716711
avg_resp_model:lang_policy3limited integration 0.507779
avg_resp_model:lang_policy3moderate or extensive 0.224106
(Intercept) .
language_profileConsistent Indigenous
language_profileMixed/Code-Switcher
avg_diff_model
avg_resp_model *
service_utilization .
Interpersonal_Index_dm ***
Institutional_Index_dm ***
age **
civic_participation_score ***
gndrwoman
urban_ruralurban ***
lvl_eduPrimary
lvl_eduSecondary **
lvl_eduHigher Education
lvl_eduNo Answer
religionMuslim ***
religionProtestant/Christian
religionOther/No Answer
religionCatholic *
lpi_dm
lang_policy3limited integration
lang_policy3moderate or extensive
ctry_Institutional_Index
ctry_Interpersonal_Index
ctry_lpi *
language_profileConsistent Indigenous:age *
language_profileMixed/Code-Switcher:age
avg_diff_model:age
avg_resp_model:age *
service_utilization:age
Interpersonal_Index_dm:age
Institutional_Index_dm:age
language_profileConsistent Indigenous:lang_policy3limited integration ***
language_profileMixed/Code-Switcher:lang_policy3limited integration
language_profileConsistent Indigenous:lang_policy3moderate or extensive
language_profileMixed/Code-Switcher:lang_policy3moderate or extensive
avg_resp_model:lang_policy3limited integration
avg_resp_model:lang_policy3moderate or extensive
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Correlation matrix not shown by default, as p = 39 > 12.
Use print(x, correlation=TRUE) or
vcov(x) if you need itanova(model_8, model_10)Data: df_mlvl
Models:
model_8: hesitancy ~ language_profile + avg_diff_model + avg_resp_model + service_utilization + Interpersonal_Index_dm + Institutional_Index_dm + civic_participation_score + gndr + age + urban_rural + lvl_edu + religion + lpi_dm + lang_policy3 + ctry_Institutional_Index + ctry_Interpersonal_Index + ctry_lpi + language_profile:lang_policy3 + avg_resp_model:lang_policy3 + (1 | ctry)
model_10: hesitancy ~ (language_profile + avg_diff_model + avg_resp_model + service_utilization + Interpersonal_Index_dm + Institutional_Index_dm) * age + civic_participation_score + gndr + urban_rural + lvl_edu + religion + lpi_dm + lang_policy3 + ctry_Institutional_Index + ctry_Interpersonal_Index + ctry_lpi + language_profile:lang_policy3 + avg_resp_model:lang_policy3 + (1 | ctry)
npar AIC BIC logLik -2*log(L) Chisq Df Pr(>Chisq)
model_8 33 28527 28792 -14230 28461
model_10 40 28529 28850 -14225 28449 11.824 7 0.1065pred_cat <- ggpredict(model_9,
terms = c("avg_resp_model", "age_group"))
p_cat <- plot(pred_cat) +
labs(title = "(a) Categorical age specification",
x = "Perceived institutional respect (standardized)",
y = "Predicted probability of hesitancy",
colour = "Age group",
fill = "Age group") +
scale_y_continuous(labels = scales::percent_format(accuracy = 1)) +
theme_minimal(base_size = 11) +
theme(legend.position = "bottom")Scale for y is already present.
Adding another scale for y, which will replace the existing scale.pred_cont <- ggpredict(model_10,
terms = c("avg_resp_model",
"age"))
p_cont <- plot(pred_cont) +
labs(title = "(b) Continuous age specification",
x = "Perceived institutional respect (standardized)",
y = "Predicted probability of hesitancy",
colour = "Age (SD units)",
fill = "Age (SD units)") +
scale_y_continuous(labels = scales::percent_format(accuracy = 1)) +
theme_minimal(base_size = 11) +
theme(legend.position = "bottom")Scale for y is already present.
Adding another scale for y, which will replace the existing scale.age_comparison_plot <- p_cat + p_cont +
plot_layout(ncol = 2) +
plot_annotation(
caption = paste(
"Note. Both specifications yield qualitatively identical patterns:",
"the protective effect of perceived institutional respect on hesitancy",
"is most pronounced among younger respondents and attenuates with age."
),
theme = theme(plot.caption = element_text(hjust = 0, size = 9,
face = "italic"))
)
# Display
print(age_comparison_plot)
isSingular(model_8)[1] FALSEisSingular(model_9)[1] FALSEmodel_performance(model_8)# Indices of model performance
AIC | AICc | BIC | R2 (cond.) | R2 (marg.) | ICC | RMSE | Sigma
-----------------------------------------------------------------------------
28527.0 | 28527.1 | 28791.8 | 0.145 | 0.060 | 0.091 | 0.467 | 1
AIC | Log_loss | Score_log
------------------------------
28527.0 | 0.626 | -Infvif(model_8) GVIF Df GVIF^(1/(2*Df))
language_profile 6.004161 2 1.565356
avg_diff_model 3.539319 1 1.881308
avg_resp_model 5.051831 1 2.247628
service_utilization 6.533634 1 2.556097
Interpersonal_Index_dm 1.089057 1 1.043579
Institutional_Index_dm 1.108358 1 1.052786
civic_participation_score 1.093454 1 1.045684
gndr 1.076388 1 1.037491
age 1.108914 1 1.053050
urban_rural 1.170289 1 1.081799
lvl_edu 1.437688 4 1.046425
religion 1.085396 4 1.010296
lpi_dm 1.079658 1 1.039066
lang_policy3 2.406060 2 1.245451
ctry_Institutional_Index 1.840938 1 1.356812
ctry_Interpersonal_Index 1.864526 1 1.365477
ctry_lpi 1.182826 1 1.087578
language_profile:lang_policy3 11.253417 4 1.353351
avg_resp_model:lang_policy3 2.780245 2 1.291281vif(model_9) GVIF Df GVIF^(1/(2*Df))
language_profile 8.378840 2 1.701359
avg_diff_model 6.032992 1 2.456215
avg_resp_model 7.151188 1 2.674171
service_utilization 11.043430 1 3.323166
Interpersonal_Index_dm 1.826258 1 1.351391
Institutional_Index_dm 1.936289 1 1.391506
age_group 8.669992 1 2.944485
civic_participation_score 1.095878 1 1.046842
gndr 1.072704 1 1.035714
urban_rural 1.171641 1 1.082424
lvl_edu 1.428441 4 1.045581
religion 1.084639 4 1.010208
lpi_dm 1.079892 1 1.039179
lang_policy3 2.405699 2 1.245404
ctry_Institutional_Index 1.840929 1 1.356808
ctry_Interpersonal_Index 1.864433 1 1.365442
ctry_lpi 1.182638 1 1.087491
language_profile:age_group 12.895274 2 1.894993
avg_diff_model:age_group 5.949469 1 2.439153
avg_resp_model:age_group 5.528351 1 2.351245
service_utilization:age_group 11.369904 1 3.371929
Interpersonal_Index_dm:age_group 1.815527 1 1.347415
Institutional_Index_dm:age_group 1.930495 1 1.389423
language_profile:lang_policy3 11.301381 4 1.354071
avg_resp_model:lang_policy3 2.782839 2 1.291582