Re: R Inla and Balancing data set
<[email protected]> Mon, 20 Apr 2026 11:10:43 -0400
| Newsgroups | gmane.comp.lang.r.general |
|---|---|
| Message-ID | <[email protected]> |
You would need to explain what you mean in more concrete terms. What, specifically, does it mean to BALANCE data" There are many ways people examine their data to test for things like validity of entries as in a column containing only positive values below 100. They may identify rows with required values that are missing and either fix them or remove them. Until we know what it is you want, we are not helpful. -----Original Message----- From: R-help <[email protected]> On Behalf Of Nasra Mapoy Sent: Sunday, April 19, 2026 4:57 PM To: [email protected] Subject: [R] R Inla and Balancing data set 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.