Re: R Inla and Balancing data set
Michael Dewey via R-help <[email protected]> Mon, 20 Apr 2026 09:39:47 +0100
| Newsgroups | gmane.comp.lang.r.general |
|---|---|
| Message-ID | <[email protected]> |
Dear Nasra
I think it is unlikely that you are going to find anyone on this list
who is willing to go through pages of your code to try to help you. You
would stand a better chance if you could produce a minimal example of
the core problem. You say you want to balance your data set so you can
give a toy example with a small data-set and as a few packages as possible.
Michael
On 19/04/2026 21:57, Nasra Mapoy wrote:
> Hi,
> My name is Nasra. I am working on my thesis using R for analysis, I have an
> objective of predicting household inequality, using R inla and machine
> learning. But my data is imbalanced (i have about 94 inadequate housing
> coded as 1) i want to balance then proceed with analysis. Unfortunately I
> got stuck each time. Please help me.
> Down is just an example of my data set
>
> head(df_analysis)# A tibble: 6 × 40
> H_REGION_NAME H_DISTRICT_NAME H_WARD_NAME H_VILLAGE_NAME H_LATITUDE
> <fct> <fct> <fct> <fct> <dbl>1
> Kusini Pemba Chake Chake Chanjaani Vijijini -5.262
> Kusini Pemba Chake Chake Chanjaani Vijijini -5.263
> Kusini Pemba Chake Chake Chanjaani Vijijini -5.264
> Kusini Pemba Chake Chake Chanjaani Vijijini -5.265
> Kusini Pemba Chake Chake Chanjaani Vijijini -5.266
> Kusini Pemba Chake Chake Chanjaani Vijijini -5.26# ℹ
> 35 more variables: H_LONGITUDE <dbl>, LOCATION <fct>, MEMBERS <dbl>,#
> TENURE <fct>, LAND <fct>, ROOF <fct>, FLOOR <fct>, WALLS <fct>,#
> ROOMS <dbl>, DRINKING_WATER <fct>, COOKING <fct>, LIGHTING <fct>,#
> TOILET <fct>, SOLID_WASTE <fct>, HEAD_SEX <fct>, HEAD_AGE <dbl>,#
> HEAD_EDUCATION <fct>, HEAD_ECON_ACTV <fct>, HEAD_MARITAL <fct>,#
> HEAD_OCCUPATION <fct>, H_Education <fct>, Durability_Floor <fct>,#
> Durability_Wall <fct>, Durability_Roof <fct>, Improved_Water <fct>,
>
> these ara some of the codes which i run
>
> library(spldv)
> library(readr)
> library(spdep)
> library(spDataLarge)
> library(haven)
> library(sp)
> library(INLA) # INLA models
> library(sf) # spatial data
> library(ggplot2) # plotting
> library(dplyr)
> #library(DMwR)
> library(smotefamily)
>
> library(sf)
> library(sp)
> library(INLA)
>
> # 1. Read shapefile and prepare coordinates
> shp_data <- st_read("housing.shp")
> shp_data_sp <- as(shp_data, "Spatial")
> coords <- st_coordinates(shp_data)
>
>
>
> # Add coordinates to dataframe
> shp_data$X <- coords[,1]
> shp_data$Y <- coords[,2]
>
> names(shp_data)[names(shp_data) == "slum"] <- "housing_inadequacy"
>
> # Convert categorical variables to factors
> shp_data$LOCATION <- as.factor(shp_data$LOCATIO)
> shp_data$HEAD_SEX <- as.factor(shp_data$HEAD_SE)
> shp_data$age_grouped <- as.factor(shp_data$ag_grpd)
> shp_data$HEAD_MARITAL<- as.factor(shp_data$HEAD_MA)
> shp_data$HEAD_EDUCATION <- as.factor(shp_data$H_Edctn)
> shp_data$HEAD_ECON_ACTV <- as.factor(shp_data$HEAD_EC)
> shp_data$hhsize <- as.factor(shp_data$hhsize)
>
> # Create numeric district id
> shp_data$id <- as.numeric(as.factor(shp_data$H_DISTR))
>
> # Response variable (binary 0/1)
> shp_data$y <- as.integer(shp_data$housing_inadequacy)
>
>
> coords <- st_coordinates(st_centroid(shp_data))
>
>
> # 2. Build mesh
> mesh <- inla.mesh.2d(
> loc = coords,
> max.edge = c(0.05, 0.2),
> cutoff = 0.01
> )
>
>
> str(coords)
>
> # 3. Define SPDE model
> spde <- inla.spde2.pcmatern(
> mesh = mesh,
> alpha = 2,
> prior.range = c(0.5, 0.01), # P(range < 0.5) = 0.01
> prior.sigma = c(1, 0.01) # P(sigma > 1) = 0.01
> )
>
> # 4. Create spatial index and projector matrix
> spatial_index <- inla.spde.make.index("spatial_field", n.spde = spde$n.spde)
>
> A <- inla.spde.make.A(mesh = mesh, loc = coords)
>
> # 5. Build stack
> stack <- inla.stack(
> data = list(y = shp_data$y),
> A = list(1, A),
> effects = list(
> data.frame(
> Intercept = 1,
> LOCATION = shp_data$LOCATION,
> age_grouped = shp_data$age_grouped,
> HEAD_SEX = shp_data$HEAD_SEX,
> HEAD_EDUCATION = shp_data$HEAD_EDUCATION,
> hhsize = shp_data$hhsize,
> HEAD_ECON_ACTV = shp_data$HEAD_ECON_ACTV,
> HEAD_MARITAL = shp_data$HEAD_MARITAL,
> id = shp_data$id
> ),
> spatial_index
> ),
> tag = "est"
> )
>
> # 6. Define formulas
>
> ## Baseline (no random effects)
> formula_baseline <- y ~ LOCATION + age_grouped + HEAD_SEX +
> HEAD_EDUCATION + hhsize + HEAD_ECON_ACTV + HEAD_MARITAL
>
> ## Unstructured (IID only)
> formula_unstructured <- y ~ LOCATION + age_grouped + HEAD_SEX +
> HEAD_EDUCATION + hhsize + HEAD_ECON_ACTV + HEAD_MARITAL +
> f(id, model = "iid")
>
> ## Structured (SPDE only)
> formula_structured <- y ~ LOCATION + age_grouped + HEAD_SEX +
> HEAD_EDUCATION + hhsize + HEAD_ECON_ACTV + HEAD_MARITAL +
> f(spatial_field, model = spde)
>
> ## Full (SPDE + IID)
> formula_full <- y ~ LOCATION + age_grouped + HEAD_SEX +
> HEAD_EDUCATION + hhsize + HEAD_ECON_ACTV + HEAD_MARITAL +
> f(id, model = "iid") +
> f(spatial_field, model = spde)
>
> # 7. Fit models
>
> result_baseline <- inla(formula_baseline, family = "binomial",
> data = inla.stack.data(stack),
> control.predictor = list(A =
> inla.stack.A(stack), compute = TRUE),
> control.compute = list(dic = TRUE, waic = TRUE))
>
> result_unstructured <- inla(formula_unstructured, family = "binomial",
> data = inla.stack.data(stack),
> control.predictor = list(A =
> inla.stack.A(stack), compute = TRUE),
> control.compute = list(dic = TRUE, waic = TRUE))
>
> result_structured <- inla(formula_structured, family = "binomial",
> data = inla.stack.data(stack),
> control.predictor = list(A =
> inla.stack.A(stack), compute = TRUE),
> control.compute = list(dic = TRUE, waic = TRUE))
>
> result_full <- inla(formula_full, family = "binomial",
> data = inla.stack.data(stack),
> control.predictor = list(A = inla.stack.A(stack),
> compute = TRUE),
> control.compute = list(dic = TRUE, waic = TRUE))
>
>
> # 8. Compare models
> model_comparison <- data.frame(
> Model = c("Baseline", "Unstructured", "Structured", "Full"),
> DIC = c(result_baseline$dic$dic,
> result_unstructured$dic$dic,
> result_structured$dic$dic,
> result_full$dic$dic),
> WAIC = c(result_baseline$waic$waic,
> result_unstructured$waic$waic,
> result_structured$waic$waic,
> result_full$waic$waic)
> )
>
> print(model_comparison)
>
>
> # Baseline model predictions
> pred_baseline <- result_baseline$summary.fitted.values
>
> # Unstructured model predictions
> pred_unstructured <- result_unstructured$summary.fitted.values
>
> # Structured model predictions
> pred_structured <- result_structured$summary.fitted.values
>
> # Full model predictions
> pred_full <- result_full$summary.fitted.values
>
> # Add posterior mean predictions
> # Extract fitted values for each model
> fv_baseline <- result_baseline$summary.fitted.values
> fv_unstructured <- result_unstructured$summary.fitted.values
> fv_structured <- result_structured$summary.fitted.values
> fv_full <- result_full$summary.fitted.values
>
> # Align with estimation stack
> idx <- inla.stack.index(stack, tag = "est")$data
>
> # Posterior mean predictions
> shp_data$pred_baseline <- fv_baseline$mean[idx]
> shp_data$pred_unstructured <- fv_unstructured$mean[idx]
> shp_data$pred_structured <- fv_structured$mean[idx]
> shp_data$pred_full <- fv_full$mean[idx]
>
> # Credible intervals
> shp_data$pred_baseline_lower <- fv_baseline$`0.025quant`[idx]
> shp_data$pred_baseline_upper <- fv_baseline$`0.975quant`[idx]
>
> shp_data$pred_unstructured_lower <- fv_unstructured$`0.025quant`[idx]
> shp_data$pred_unstructured_upper <- fv_unstructured$`0.975quant`[idx]
>
> shp_data$pred_structured_lower <- fv_structured$`0.025quant`[idx]
> shp_data$pred_structured_upper <- fv_structured$`0.975quant`[idx]
>
> shp_data$pred_full_lower <- fv_full$`0.025quant`[idx]
> shp_data$pred_full_upper <- fv_full$`0.975quant`[idx]
>
>
> ##Plotting
> library(sf)
> library(ggplot2)
> library(patchwork) # for arranging plots
>
> tanzania_boundaries <- st_read("C:/Users/Malilo/Desktop/New
> folder/PHD_2026/TPHC_10/tza_admbnda_adm1_20181019/tza_admbnda_adm1_20181019.shp")
>
>
> # Create centroids for region labels
> library(sf)
> library(ggplot2)
> library(viridis)
> library(ggspatial) # for north arrow and scale bar
>
>
> # Create centroids for labels
> tanzania_centroids <- st_centroid(tanzania_boundaries)
>
> # Convert back to sf object
> shp_data_sf <- st_as_sf(shp_data)
>
> # Baseline model
> p_baseline <- ggplot() +
> geom_sf(data = tanzania_boundaries, fill = NA, color = "black",
> linewidth = 0.6) +
> geom_sf(data = shp_data_sf, aes(color = pred_baseline)) +
> geom_sf_text(data = tanzania_centroids, aes(label = ADM1_EN),
> size = 2, color = "black", check_overlap = TRUE) +
> scale_color_viridis_c(option = "magma") +
> theme_minimal() +
> theme(
> panel.grid.major = element_blank(),
> panel.grid.minor = element_blank(),
> axis.text = element_blank(),
> axis.ticks = element_blank(),
> axis.title = element_blank(),
> panel.border = element_rect(color = "black", fill = NA, linewidth = 1)
> ) +
> labs(title = "Housing Inadequacy (Baseline Model)", color = "Probability")
>
> print(p_baseline)
>
>
> # Unstructured model
> p_unstructured <- ggplot() +
> geom_sf(data = tanzania_boundaries, fill = NA, color = "black",
> linewidth = 0.6) +
> geom_sf(data = shp_data_sf, aes(color = pred_unstructured)) +
> geom_sf_text(data = tanzania_centroids, aes(label = ADM1_EN),
> size = 2, color = "black", check_overlap = TRUE) +
> scale_color_viridis_c(option = "magma") +
> theme_minimal() +
> theme(
> panel.grid.major = element_blank(),
> panel.grid.minor = element_blank(),
> axis.text = element_blank(),
> axis.ticks = element_blank(),
> axis.title = element_blank(),
> panel.border = element_rect(color = "black", fill = NA, linewidth = 1)
> ) +
> labs(title = "Housing Inadequacy (Unstructured Model)", color = "Probability")
>
> print(p_unstructured)
>
>
> # Structured model
> p_structured <- ggplot() +
> geom_sf(data = tanzania_boundaries, fill = NA, color = "black",
> linewidth = 0.6) +
> geom_sf(data = shp_data_sf, aes(color = pred_structured)) +
> geom_sf_text(data = tanzania_centroids, aes(label = ADM1_EN),
> size = 2, color = "black", check_overlap = TRUE) +
> scale_color_viridis_c(option = "magma") +
> theme_minimal() +
> theme(
> panel.grid.major = element_blank(),
> panel.grid.minor = element_blank(),
> axis.text = element_blank(),
> axis.ticks = element_blank(),
> axis.title = element_blank(),
> panel.border = element_rect(color = "black", fill = NA, linewidth = 1)
> ) +
> labs(title = "Housing Inadequacy (Structured Model)", color = "Probability")
>
> print(p_structured)
>
>
> # Full model
> p_full <- ggplot() +
> geom_sf(data = tanzania_boundaries, fill = NA, color = "black",
> linewidth = 0.6) +
> geom_sf(data = shp_data_sf, aes(color = pred_full)) +
> geom_sf_text(data = tanzania_centroids, aes(label = ADM1_EN),
> size = 2, color = "black", check_overlap = TRUE) +
> scale_color_viridis_c(option = "magma") +
> theme_minimal() +
> theme(
> panel.grid.major = element_blank(),
> panel.grid.minor = element_blank(),
> axis.text = element_blank(),
> axis.ticks = element_blank(),
> axis.title = element_blank(),
> panel.border = element_rect(color = "black", fill = NA, linewidth = 1)
> ) +
> labs(title = "Housing Inadequacy (Full Model)", color = "Probability")
>
> print(p_full)
>
>
> # Define mapping of coded names to descriptive labels
> label_map <- c(
> "(Intercept)" = "(Intercept)",
> "LOCATION2" = "Location: Urban",
> "HEAD_SEX2" = "Head sex: Female",
> "age_grouped2"= "Age group: 15–24",
> "age_grouped3"= "Age group: 25–44",
> "age_grouped4"= "Age group: 45–64",
> "age_grouped5"= "Age group: 65+",
> "hhsize2" = "Household size: 2 - 3 persons",
> "hhsize3" = "Household size: 4 - 5 persons",
> "hhsize4" = "Household size: 6 persons or above",
> "HEAD_MARITAL2"= "Marital Status: Married",
> "HEAD_MARITAL3"= "Marital Status: Living together",
> "HEAD_MARITAL4"= "Marital Status: Divorced",
> "HEAD_MARITAL5"= "Marital Status: Separated",
> "HEAD_MARITAL6"= "Marital Status: Widowed",
> "HEAD_EDUCATION2"= "Education: Secondary",
> "HEAD_EDUCATION3"= "Education: Tertiary",
> "HEAD_ECON_ACTV2"= "Self-employed (No workers) ",
> "HEAD_ECON_ACTV3"= "Self-employed (With workers) ",
> "HEAD_ECON_ACTV4"= "Unskilled labour",
> "HEAD_ECON_ACTV5"= "Unrecognized employees"
> )
>
>
>
> # Function to relabel rownames
> relabel_fixed <- function(fixed_table, map){
> fixed_table$Predictor <- map[rownames(fixed_table)]
> rownames(fixed_table) <- NULL
> return(fixed_table)
> }
>
> fixed_baseline_lbl <- relabel_fixed(result_baseline$summary.fixed,
> label_map)
> fixed_unstructured_lbl <-
> relabel_fixed(result_unstructured$summary.fixed, label_map)
> fixed_structured_lbl <-
> relabel_fixed(result_structured$summary.fixed, label_map)
> fixed_full_lbl <- relabel_fixed(result_full$summary.fixed, label_map)
>
>
> library(openxlsx)
>
> wb <- createWorkbook()
> addWorksheet(wb, "Baseline")
> addWorksheet(wb, "Unstructured")
> addWorksheet(wb, "Structured")
> addWorksheet(wb, "Full")
>
> writeData(wb, "Baseline", fixed_baseline_lbl)
> writeData(wb, "Unstructured", fixed_unstructured_lbl)
> writeData(wb, "Structured", fixed_structured_lbl)
> writeData(wb, "Full", fixed_full_lbl)
>
> saveWorkbook(wb, "Fixed_Effects_Models_Labeled.xlsx", overwrite = TRUE)
>
>
>
> library(ggplot2)
>
> # Relabel function
> relabel_for_plot <- function(fixed_table, map){
> df <- fixed_table
> df$Predictor <- map[rownames(df)]
> rownames(df) <- NULL
> return(df)
> }
>
> # Relabel each model’s fixed effects
> fixed_baseline_lbl <-
> relabel_for_plot(result_baseline$summary.fixed, label_map)
> fixed_unstructured_lbl <-
> relabel_for_plot(result_unstructured$summary.fixed, label_map)
> fixed_structured_lbl <-
> relabel_for_plot(result_structured$summary.fixed, label_map)
> fixed_full_lbl <- relabel_for_plot(result_full$summary.fixed, label_map)
>
> # Plot Baseline
> ggplot(fixed_baseline_lbl, aes(y = Predictor, x = mean)) +
> geom_point(color = "blue") +
> geom_errorbarh(aes(xmin = `0.025quant`, xmax = `0.975quant`), height = 0.2) +
> theme_minimal() +
> labs(title = "Posterior Estimates with 95% Credible Intervals
> (Baseline Model)",
> y = "Predictors", x = "Posterior Mean")
>
> # Plot Unstructured
> ggplot(fixed_unstructured_lbl, aes(y = Predictor, x = mean)) +
> geom_point(color = "blue") +
> geom_errorbarh(aes(xmin = `0.025quant`, xmax = `0.975quant`), height = 0.2) +
> theme_minimal() +
> labs(title = "Posterior Estimates with 95% Credible Intervals
> (Unstructured Model)",
> y = "Predictors", x = "Posterior Mean")
>
> # Plot Structured
> ggplot(fixed_structured_lbl, aes(y = Predictor, x = mean)) +
> geom_point(color = "blue") +
> geom_errorbarh(aes(xmin = `0.025quant`, xmax = `0.975quant`), height = 0.2) +
> theme_minimal() +
> labs(title = "Posterior Estimates with 95% Credible Intervals
> (Structured Model)",
> y = "Predictors", x = "Posterior Mean")
>
> # Plot Full
> ggplot(fixed_full_lbl, aes(y = Predictor, x = mean)) +
> geom_point(color = "blue") +
> geom_errorbarh(aes(xmin = `0.025quant`, xmax = `0.975quant`), height = 0.2) +
> theme_minimal() +
> labs(title = "Posterior Estimates with 95% Credible Intervals (Full Model)",
> y = "Predictors", x = "Posterior Mean")
>
>
> library(pROC)
>
> # Example for Full model
> roc_full <- roc(response = shp_data_sf$y, predictor = shp_data_sf$pred_full)
>
> plot(roc_full, col = "blue", main = "ROC Curve - Full Model")
> auc(roc_full) # Print AUC value
>
> # Compare across models
> roc_baseline <- roc(response = shp_data_sf$y, predictor =
> shp_data_sf$pred_baseline)
> roc_unstructured <- roc(response = shp_data_sf$y, predictor =
> shp_data_sf$pred_unstructured)
> roc_structured <- roc(response = shp_data_sf$y, predictor =
> shp_data_sf$pred_structured)
>
>
> plot(roc_baseline, col = "red", main = "ROC Curves Comparison")
> plot(roc_unstructured, col = "green", add = TRUE)
> plot(roc_structured, col = "purple", add = TRUE)
> plot(roc_full, col = "blue", add = TRUE)
> legend("bottomright", legend = c("Baseline","Unstructured","Structured","Full"),
> col = c("red","green","purple","blue"), lwd = 2)
>
>
>
> ##Baseline Model
> # Fixed effects
> fixed_baseline_lbl <- relabel_for_plot(result_baseline$summary.fixed, label_map)
>
> ggplot(fixed_baseline_lbl, aes(y = Predictor, x = mean)) +
> geom_point(color = "blue") +
> geom_errorbarh(aes(xmin = `0.025quant`, xmax = `0.975quant`), height = 0.2) +
> labs(title = "Posterior Estimates (Baseline Model)", y =
> "Predictors", x = "Posterior Mean") +
> theme_minimal()
>
> # Predicted probabilities (per observation)
> pred_df_baseline <- data.frame(ID = 1:nrow(shp_data),
> mean = shp_data$pred_baseline)
> ggplot(pred_df_baseline, aes(x = ID, y = mean)) +
> geom_point(color = "blue") +
> labs(title = "Predicted Probabilities (Baseline Model)",
> x = "Observation Index", y = "Posterior Mean") +
> theme_minimal()
>
> # Histogram of predictions
> ggplot(pred_df_baseline, aes(x = mean)) +
> geom_histogram(binwidth = 0.05, fill = "skyblue", color = "black") +
> labs(title = "Distribution of Predicted Probabilities (Baseline Model)",
> x = "Predicted Probability", y = "Count") +
> theme_minimal()
>
> # ROC curve
> roc_baseline <- roc(response = shp_data$y, predictor = shp_data$pred_baseline)
> plot(roc_baseline, col = "red", main = "ROC Curve (Baseline)")
> auc(roc_baseline)
>
>
> ##Unstructured Model
> # Fixed effects
> fixed_unstructured_lbl <-
> relabel_for_plot(result_unstructured$summary.fixed, label_map)
>
> ggplot(fixed_unstructured_lbl, aes(y = Predictor, x = mean)) +
> geom_point(color = "blue") +
> geom_errorbarh(aes(xmin = `0.025quant`, xmax = `0.975quant`), height = 0.2) +
> labs(title = "Posterior Estimates (Unstructured Model)", y =
> "Predictors", x = "Posterior Mean") +
> theme_minimal()
>
> # Predicted probabilities (per observation)
> pred_df_unstructured <- data.frame(ID = 1:nrow(shp_data),
> mean = shp_data$pred_unstructured)
> ggplot(pred_df_unstructured, aes(x = ID, y = mean)) +
> geom_point(color = "red") +
> labs(title = "Predicted Probabilities (Unstructured Model)",
> x = "Observation Index", y = "Posterior Mean") +
> theme_minimal()
>
> # Histogram
> ggplot(pred_df_unstructured, aes(x = mean)) +
> geom_histogram(binwidth = 0.05, fill = "salmon", color = "black") +
> labs(title = "Distribution of Predicted Probabilities (Unstructured Model)",
> x = "Predicted Probability", y = "Count") +
> theme_minimal()
>
> # Random effects (IID)
> rand_df_unstructured <- as.data.frame(result_unstructured$summary.random$id)
> rand_df_unstructured$Index <- 1:nrow(rand_df_unstructured)
>
> ggplot(rand_df_unstructured, aes(x = Index, y = mean)) +
> geom_point(color = "blue") +
> labs(title = "Unstructured Random Effect (IID)",
> x = "District Index", y = "Posterior Mean") +
> theme_minimal()
>
>
> # Extract IID random effects
> # Convert shapefile data to sf object
> shp_data_sf <- st_as_sf(shp_data)
>
> # Map predicted probabilities from the Unstructured model
> ggplot() +
> geom_sf(data = tanzania_boundaries, fill = NA, color = "black",
> linewidth = 0.6) +
> geom_sf(data = shp_data_sf, aes(color = pred_unstructured)) +
> geom_sf_text(data = tanzania_centroids, aes(label = ADM1_EN),
> size = 2, color = "black", check_overlap = TRUE) +
> scale_color_viridis_c(option = "magma") +
> theme_minimal() +
> theme(
> panel.grid.major = element_blank(),
> panel.grid.minor = element_blank(),
> axis.text = element_blank(),
> axis.ticks = element_blank(),
> axis.title = element_blank(),
> panel.border = element_rect(color = "black", fill = NA, linewidth = 1)
> ) +
> labs(title = "Distribution of Predicted Probabilities (Unstructured Model)",
> color = "Probability")
>
> # ROC curve
> roc_unstructured <- roc(response = shp_data$y, predictor =
> shp_data$pred_unstructured)
> plot(roc_unstructured, col = "green", main = "ROC Curve (Unstructured)")
> auc(roc_unstructured)
>
> hyper_unstructured <- result_unstructured$summary.hyperpar
> print(as.data.frame(hyper_unstructured))
>
>
>
>
> ##Structured Model
> # Fixed effects
> fixed_structured_lbl <-
> relabel_for_plot(result_structured$summary.fixed, label_map)
>
> ggplot(fixed_structured_lbl, aes(y = Predictor, x = mean)) +
> geom_point(color = "blue") +
> geom_errorbarh(aes(xmin = `0.025quant`, xmax = `0.975quant`), height = 0.2) +
> labs(title = "Posterior Estimates (Structured Model)", y =
> "Predictors", x = "Posterior Mean") +
> theme_minimal()
>
> # Predicted probabilities (per observation)
> pred_df_structured <- data.frame(ID = 1:nrow(shp_data),
> mean = shp_data$pred_structured)
> ggplot(pred_df_structured, aes(x = ID, y = mean)) +
> geom_point(color = "green") +
> labs(title = "Predicted Probabilities (Structured Model)",
> x = "Observation Index", y = "Posterior Mean") +
> theme_minimal()
>
> # Histogram
> ggplot(pred_df_structured, aes(x = mean)) +
> geom_histogram(binwidth = 0.05, fill = "lightgreen", color = "black") +
> labs(title = "Distribution of Predicted Probabilities (Structured Model)",
> x = "Predicted Probability", y = "Count") +
> theme_minimal()
>
> # Random effects (SPDE field)
> rand_df_structured <-
> as.data.frame(result_structured$summary.random$spatial_field)
> rand_df_structured$Index <- 1:nrow(rand_df_structured)
>
> ggplot(rand_df_structured, aes(x = Index, y = mean)) +
> geom_point(color = "darkgreen") +
> labs(title = "Structured Random Effect (SPDE)",
> x = "Spatial Index", y = "Posterior Mean") +
> theme_minimal()
> # Map predicted probabilities (Structured)
> ggplot() +
> geom_sf(data = tanzania_boundaries, fill = NA, color = "black",
> linewidth = 0.6) +
> geom_sf(data = shp_data_sf, aes(color = pred_structured)) +
> geom_sf_text(data = tanzania_centroids, aes(label = ADM1_EN),
> size = 2, color = "black", check_overlap = TRUE) +
> scale_color_viridis_c(option = "magma") +
> labs(title = "Predicted Probabilities (Structured Model)", color =
> "Probability") +
> theme_minimal()
>
>
> # ROC curve
> roc_structured <- roc(response = shp_data$y, predictor =
> shp_data$pred_structured)
> plot(roc_structured, col = "purple", main = "ROC Curve (Structured)")
> auc(roc_structured)
>
> # Extract hyperparameters for Structured model
> hyper_structured <- result_structured$summary.hyperpar
>
> # Convert to data frame
> hyper_structured_df <- as.data.frame(hyper_structured)
>
> # Add descriptive row names
> rownames(hyper_structured_df) <- c("Theta1 for spatial", "Theta2 for spatial")
>
> # Print table
> print(hyper_structured_df)
>
> # Optional: export to Excel
> library(openxlsx)
> wb <- createWorkbook()
> addWorksheet(wb, "Structured Hyperparameters")
> writeData(wb, "Structured Hyperparameters", hyper_structured_df,
> rowNames = TRUE)
> saveWorkbook(wb, "Structured_Hyperparameters.xlsx", overwrite = TRUE)
>
>
>
>
>
> ##Full Model
> # Fixed effects
> fixed_full_lbl <- relabel_for_plot(result_full$summary.fixed, label_map)
>
> ggplot(fixed_full_lbl, aes(y = Predictor, x = mean)) +
> geom_point(color = "blue") +
> geom_errorbarh(aes(xmin = `0.025quant`, xmax = `0.975quant`), height = 0.2) +
> labs(title = "Posterior Estimates (Full Model)", y = "Predictors", x
> = "Posterior Mean") +
> theme_minimal()
>
> # Predicted probabilities (per observation)
> pred_df_full <- data.frame(ID = 1:nrow(shp_data),
> mean = shp_data$pred_full)
> ggplot(pred_df_full, aes(x = ID, y = mean)) +
> geom_point(color = "black") +
> labs(title = "Predicted Probabilities (Full Model)",
> x = "Observation Index", y = "Posterior Mean") +
> theme_minimal()
>
> # Histogram
> ggplot(pred_df_full, aes(x = mean)) +
> geom_histogram(binwidth = 0.05, fill = "gray", color = "black") +
> labs(title = "Distribution of Predicted Probabilities (Full Model)",
> x = "Predicted Probability", y = "Count") +
> theme_minimal()
>
> # Random effects (IID)
> rand_df_full_iid <- as.data.frame(result_full$summary.random$id)
> rand_df_full_iid$Index <- 1:nrow(rand_df_full_iid)
>
> ggplot(rand_df_full_iid, aes(x = Index, y = mean)) +
> geom_point(color = "purple") +
> labs(title = "Full Model Random Effect (IID)",
> x = "District Index", y = "Posterior Mean") +
> theme_minimal()
>
> # Random effects (SPDE)
> rand_df_full_spatial <- as.data.frame(result_full$summary.random$spatial_field)
> rand_df_full_spatial$Index <- 1:nrow(rand_df_full_spatial)
>
> ggplot(rand_df_full_spatial, aes(x = Index, y = mean)) +
> geom_point(color = "orange") +
> labs(title = "Full Model Random Effect (SPDE)",
> x = "Spatial Index", y = "Posterior Mean") +
> theme_minimal()
>
> # Map predicted probabilities (Full)
> ggplot() +
> geom_sf(data = tanzania_boundaries, fill = NA, color = "black",
> linewidth = 0.6) +
> geom_sf(data = shp_data_sf, aes(color = pred_full)) +
> geom_sf_text(data = tanzania_centroids, aes(label = ADM1_EN),
> size = 2, color = "black", check_overlap = TRUE) +
> scale_color_viridis_c(option = "magma") +
> labs(title = "Predicted Probabilities (Full Model)", color = "Probability") +
> theme_minimal()
>
>
> # ROC curve
> roc_full <- roc(response = shp_data$y, predictor = shp_data$pred_full)
> plot(roc_full, col = "blue", main = "ROC Curve (Full)")
> auc(roc_full)
>
> hyper_full <- result_full$summary.hyperpar
> print(as.data.frame(hyper_full))
>
>
> library(pROC)
> library(yardstick)
> library(dplyr)
>
> # Function to compute metrics for any INLA model
> evaluate_inla <- function(model, stack, df_inla, model_name) {
>
> # Prediction indices
> index_pred <- inla.stack.index(stack, "pred")$data
>
> # True labels
> true_labels <- df_inla$y[index_pred]
>
> # Posterior predicted probabilities
> pred_probs <- model$summary.fitted.values$mean[index_pred]
>
> # ROC + AUC
> roc_obj <- roc(true_labels, pred_probs)
> auc_val <- auc(roc_obj)
>
> # Optimal threshold (Youden’s J)
> opt_thresh <- coords(roc_obj, "best", ret = "threshold")
>
> # Classify
> pred_class <- ifelse(pred_probs > opt_thresh, 1, 0)
>
> # Evaluation dataframe
> df_eval <- data.frame(
> truth = factor(true_labels, levels = c(0,1)),
> pred_class = factor(pred_class, levels = c(0,1)),
> pred_prob = pred_probs
> )
>
> # Metrics
> acc <- accuracy(df_eval, truth, pred_class)$.estimate
> prec <- precision(df_eval, truth, pred_class)$.estimate
> rec <- recall(df_eval, truth, pred_class)$.estimate
> f1 <- f_meas(df_eval, truth, pred_class, beta = 1)$.estimate
>
> # Collect results
> data.frame(
> Model = model_name,
> Threshold = round(opt_thresh,3),
> Accuracy = round(acc,3),
> Precision = round(prec,3),
> Recall = round(rec,3),
> F1 = round(f1,3),
> AUC = round(auc_val,3)
> )
> }
>
> results_baseline <- evaluate_inla(result_baseline, stack,
> shp_data_sf, "INLA Baseline")
> results_unstructured<- evaluate_inla(result_unstructured, stack,
> shp_data_sf, "INLA Unstructured")
> results_structured <- evaluate_inla(result_spatial, stack,
> shp_data_sf, "INLA Structured")
> results_full <- evaluate_inla(result_full, stack, shp_data_sf,
> "INLA Full")
>
> # Combine into one table
> results_all <- bind_rows(results_baseline, results_unstructured,
> results_structured, results_full)
>
> print(results_all)
>
>
>
>
>
>
> #Machine Learning
>
> library(spldv)
> library(readr)
> library(spdep)
> library(spDataLarge)
> library(haven)
> library(sp)
> library(INLA) # INLA models
> library(sf) # spatial data
> library(ggplot2) # plotting
> library(dplyr)
> library(DMwR)
> library(smotefamily)
> library(caret)
>
> # 1. Read shapefile and prepare coordinates
> shp_data <- st_read("housing.shp")
>
> # Convert categorical variables to factors
> shp_data$LOCATION <- as.factor(shp_data$LOCATIO)
> shp_data$HEAD_SEX <- as.factor(shp_data$HEAD_SE)
> shp_data$age_grouped <- as.factor(shp_data$ag_grpd)
> shp_data$HEAD_MARITAL<- as.factor(shp_data$HEAD_MA)
> shp_data$HEAD_EDUCATION <- as.factor(shp_data$H_Edctn)
> shp_data$HEAD_ECON_ACTV <- as.factor(shp_data$HEAD_EC)
> shp_data$hhsize <- as.factor(shp_data$hhsize)
>
>
> names(shp_data)
>
>
> table(shp_data$slum)
> table(shp_data$slum)
>
>
>
> #install.packages("remotes")
> #remotes::install_version("DMwR", version = "0.4.1", repos =
> "http://cran.us.r-project.org")
>
> set.seed(123)
>
> # 1. Drop geometry column
> shp_data_no_geom <- shp_data %>% select(-geometry)
>
> # 2. Train/test split (70/30)
> trainIndex <- createDataPartition(shp_data_no_geom$slum, p = 0.7, list = FALSE)
> train <- shp_data_no_geom[trainIndex, ]
> test <- shp_data_no_geom[-trainIndex, ]
>
> # 3. Apply SMOTE only to training set (formula interface)
> # Auto-tune perc.over and perc.under for perfect balance
> minority_count <- sum(train$slum == 0)
> majority_count <- sum(train$slum == 1)
>
> perc.over <- ((majority_count - minority_count) / minority_count) * 100
> perc.under <- 100
>
> train_df <- as.data.frame(train)
> train_df$slum <- factor(train_df$slum, levels = c(0,1))
>
> # Convert all character columns to factors
> char_cols <- sapply(train_df, is.character)
> train_df[char_cols] <- lapply(train_df[char_cols], factor)
>
>
> list_cols <- sapply(train_df, is.list)
> train_df <- train_df[ , !list_cols]
>
>
> table(train_bal$slum)
>
> library(sf)
> library(dplyr)
>
> # Ensure slum is numeric in train_bal
> train_bal$slum <- as.numeric(as.character(train_bal$slum))
>
> # Join balanced attributes back to geometry
> train_bal_sf <- shp_data %>%
> inner_join(train_bal, by = c("LOCATION","age_grouped","HEAD_SEX",
> "HEAD_EDUCATION","hhsize","HEAD_ECON_ACTV",
> "HEAD_MARITAL","slum"))
>
> # Save as shapefile
> st_write(train_bal_sf, "train_balanced.shp", delete_layer = TRUE)
>
> # Load the shapefile you just wrote
> train_bal_sf <- st_read("train_balanced.shp")
>
> # Inspect column names
> names(train_bal_sf)
>
> # Check balance
> table(train_bal_sf$slum)
>
>
> names(train_balanced.shp)
>
>
> ###Modelling
>
> library(caret)
> library(randomForest)
> library(xgboost)
> library(pROC)
>
> set.seed(123)
>
> # Define baseline formula
> formula_baseline <- slum ~ LOCATION + age_grouped + HEAD_SEX +
> HEAD_EDUCATION + hhsize + HEAD_ECON_ACTV + HEAD_MARITAL
>
> # -------------------------------
> # 1. Logistic Regression (GLM)
> # -------------------------------
> logit_model <- glm(formula_baseline, data = train_bal, family = binomial)
>
> logit_pred <- predict(logit_model, newdata = test, type = "response")
> logit_class <- ifelse(logit_pred > 0.5, 1, 0)
>
> confusionMatrix(factor(logit_class), factor(test$slum))
> roc(test$slum, logit_pred)
>
> # -------------------------------
> # 2. Random Forest
> # -------------------------------
> rf_model <- randomForest(formula_baseline, data = train_bal,
> ntree = 500, mtry = 3, importance = TRUE)
>
> rf_pred <- predict(rf_model, newdata = test, type = "prob")[,2]
> rf_class <- ifelse(rf_pred > 0.5, 1, 0)
>
> confusionMatrix(factor(rf_class), factor(test$slum))
> roc(test$slum, rf_pred)
>
> # -------------------------------
> # 3. XGBoost
> # -------------------------------
> # For XGBoost we need dummy encoding
> train_X <- model.matrix(formula_baseline, data = train_bal)[, -1]
> train_y <- as.numeric(as.character(train_bal$slum))
> test_X <- model.matrix(formula_baseline, data = test)[, -1]
> test_y <- as.numeric(as.character(test$slum))
>
> xgb_train <- xgb.DMatrix(data = train_X, label = train_y)
> xgb_test <- xgb.DMatrix(data = test_X, label = test_y)
>
> xgb_model <- xgboost(data = xgb_train, nrounds = 200,
> objective = "binary:logistic", eval_metric = "auc")
>
> xgb_pred <- predict(xgb_model, newdata = xgb_test)
> xgb_class <- ifelse(xgb_pred > 0.5, 1, 0)
>
> confusionMatrix(factor(xgb_class), factor(test_y))
> roc(test_y, xgb_pred)
>
>
> Structure
>
> str(df_analysis)tibble [1,417,184 × 40] (S3: tbl_df/tbl/data.frame)
> $ H_REGION_NAME : Factor w/ 32 levels "Kusini Pemba",..: 1 1 1 1
> 1 1 1 1 1 1 ...
> $ H_DISTRICT_NAME : Factor w/ 151 levels "Chake Chake",..: 1 1 1 1
> 1 1 1 1 1 1 ...
> $ H_WARD_NAME : Factor w/ 3995 levels "Chanjaani","Shungi",..:
> 1 1 1 1 1 1 1 1 1 1 ...
> $ H_VILLAGE_NAME : Factor w/ 15570 levels "Vijijini","Muembe
> Tongwa",..: 1 1 1 1 1 1 1 1 1 1 ...
> $ H_LATITUDE : num [1:1417184] -5.26 -5.26 -5.26 -5.26 -5.26 ...
> $ H_LONGITUDE : num [1:1417184] 39.8 39.8 39.8 39.8 39.8 ...
> $ LOCATION : Factor w/ 2 levels "Urban","Rural": 1 1 1 1 1 1
> 1 1 1 1 ...
> $ MEMBERS : num [1:1417184] 10 3 5 12 1 6 7 4 8 10 ...
> $ TENURE : Factor w/ 7 levels "Owned by household",..: 1 2
> 1 1 2 2 1 2 1 1 ...
> $ LAND : Factor w/ 8 levels "Title deed","No legal
> right",..: 1 NA 2 3 NA NA 1 NA 3 2 ...
> $ ROOF : Factor w/ 8 levels "Iron sheets",..: 1 1 1 1 2
> 1 1 1 1 1 ...
> $ FLOOR : Factor w/ 10 levels "Cement","Earth/Sand",..: 1
> 1 1 1 2 1 1 1 1 1 ...
> $ WALLS : Factor w/ 10 levels "Cement bricks/Rock
> bricks",..: 1 1 1 1 1 1 1 2 1 1 ...
> $ ROOMS : num [1:1417184] 2 3 3 5 2 1 2 3 4 6 ...
> $ DRINKING_WATER : Factor w/ 14 levels "Public tap/standpipe",..:
> 1 1 2 3 1 3 2 2 3 4 ...
> $ COOKING : Factor w/ 14 levels "Firewood","Charcoal",..: 1
> 2 1 1 2 1 1 1 2 1 ...
> $ LIGHTING : Factor w/ 13 levels "Kerosene (Wick lamps)",..:
> 1 2 3 2 3 2 2 1 2 2 ...
> $ TOILET : Factor w/ 11 levels "Flush/pour flush to
> covered pit",..: 1 1 1 1 2 1 1 3 1 1 ...
> $ SOLID_WASTE : Factor w/ 9 levels "Irregularly collected",..:
> 1 2 3 3 3 3 2 4 4 4 ...
> $ HEAD_SEX : Factor w/ 2 levels "Female Headed",..: 1 2 2 2
> 2 2 2 1 2 2 ...
> $ HEAD_AGE : num [1:1417184] 44 24 34 44 24 38 47 41 64 44 ...
> $ HEAD_EDUCATION : Factor w/ 29 levels "14","10","7",..: NA 1 2 2
> 2 3 2 2 2 4 ...
> $ HEAD_ECON_ACTV : Factor w/ 5 levels "Paid employees",..: NA NA
> NA NA NA NA NA NA NA NA ...
> $ HEAD_MARITAL : Factor w/ 7 levels "Divorced","Married",..: 1 2
> 2 2 3 2 2 2 2 2 ...
> $ HEAD_OCCUPATION : Factor w/ 9 levels "Technicians and associate
> professionals.",..: 1 1 2 3 NA NA 4 NA 2 3 ...
> $ H_Education : Factor w/ 4 levels "Secondary","Primary",..: NA
> 1 1 1 1 2 1 1 1 2 ...
> $ Durability_Floor : Factor w/ 2 levels "Durable","Not durable": 1 1
> 1 1 2 1 1 1 1 1 ...
> $ Durability_Wall : Factor w/ 2 levels "Durable","Not durable": 1 1
> 1 1 1 1 1 2 1 1 ...
> $ Durability_Roof : Factor w/ 2 levels "Durable","Not durable": 1 1
> 1 1 2 1 1 1 1 1 ...
> $ Improved_Water : Factor w/ 2 levels "Improved","Unimproved": 1 1
> 1 1 1 1 1 1 1 2 ...
> $ Toilet_Facility : Factor w/ 2 levels "Unimproved","Improved": 1 1
> 1 1 2 1 1 2 1 1 ...
> $ Secure_Tenure : Factor w/ 2 levels "Not secure","Secure": 1 NA
> 2 1 NA NA 1 NA 1 2 ...
> $ age_grouped : Factor w/ 5 levels "25-44","15-24",..: 1 2 1 1
> 2 1 3 1 3 1 ...
> $ Bed_Rooms : num [1:1417184] 0 1 1 1 0 0 0 1 1 1 ...
> $ hhsize : Factor w/ 4 levels "6 and above persons",..: 1
> 2 3 1 4 1 1 3 1 1 ...
> $ Living_Area : Factor w/ 278 levels
> "5","1","1.66666662693024",..: 1 2 3 4 5 6 7 8 9 3 ...
> $ Living_Area1 : Factor w/ 2 levels "Overcrowding",..: 1 2 2 2 2
> 1 1 2 2 2 ...
> $ Durable_Materials : Factor w/ 2 levels "Durable","Not Durable": 1 1
> 1 1 2 1 1 2 1 1 ...
> $ housing_inadequacy: Factor w/ 2 levels "Inadequate","Adequate": 1 1
> 1 1 1 1 1 1 1 1 ...
> $ number_slum1 : num [1:1417184] 3 1 1 2 1 2 3 1 2 2 ...
>
>>
> sample data set
> source("sample_dataset.R")
>
>
> I am looking forward to get your feedback
> Regards
> Nasra
>
> [[alternative HTML version deleted]]
>
> ______________________________________________
> [email protected] mailing list -- To UNSUBSCRIBE and more, see
> https://stat.ethz.ch/mailman/listinfo/r-help
> PLEASE do read the posting guide https://www.R-project.org/posting-guide.html
> and provide commented, minimal, self-contained, reproducible code.
--
Michael Dewey