Understanding Methane Emissions Through Heat Map and Comparison of Many Models

It is a little-known fact that methane is more dangerous than carbon dioxide when it comes to global warming; it would harm the atmosphere about 80 times more than CO2. In this article, we will take a look at this situation across regions.

First, we examine the methane rates with a heatmap.

library(tidyverse)
library(openxlsx)
library(WDI)
library(ragg)

#Dataset for methane emissions
df_methane <- 
  WDI(indicator = "EN.ATM.METH.KT.CE", 
      extra = TRUE) %>% 
  as_tibble() %>% 
  janitor::clean_names() %>% 
  #removing labels attribute for fitting process
  crosstable::remove_labels() %>% 
  drop_na() %>% 
  rename(methane = en_atm_meth_kt_ce) 

#Dataset for heatmap  
df_heatmap <- 
  df_methane %>% 
  filter(region != "Aggregates") %>% 
  group_by(region, year) %>% 
  summarise(methane_agg = mean(methane)) %>% 
  select(region, year, methane_agg)

#Plotting heatmap for methane emissions based on region
df_heatmap %>%
  ggplot(aes(
    x = year,
    y = reorder(region, methane_agg),
    fill = methane_agg
  )) +
  geom_tile(colour = "white", linewidth = 1) +
  scale_fill_gradient(low = "steelblue",
                      high = "orange",
                      labels = scales::label_number(scale_cut = scales::cut_short_scale()),
                      name = "Methane") +
  #removing the padding between heatmap tiles and the axis
  scale_x_continuous(expand = c(0, 0), position = "top") +
  scale_y_discrete(expand = c(0, 0)) +
  labs(title = "Methane Emissions (kt of CO2 equivalent) based on Region") +
  theme(text = element_text(size = 15, family = "Bricolage Grotesque"),
        axis.title = element_blank(),
        legend.title = element_text(size = 12, face = "bold"),
        axis.ticks = element_blank(),
        plot.title = element_text(size = 15, 
                                  face = "bold",
                                  hjust = 0.5))

< svg xmlns="http://www.w3.org/2000/svg" width="12" height="12" fill="none" viewBox="0 0 12 12"> < path fill="#fff" d="M2 0a2 2 0 0 0-2 2v2h1.5V2a.5.5 0 0 1 .5-.5h2V0H2Zm2 10.5H2a.5.5 0 0 1-.5-.5V8H0v2a2 2 0 0 0 2 2h2v-1.5ZM8 12v-1.5h2a.5.5 0 0 0 .5-.5V8H12v2a2 2 0 0 1-2 2H8Zm2-12a2 2 0 0 1 2 2v2h-1.5V2a.5.5 0 0 0-.5-.5H8V0h2Z" />

< svg xmlns="http://www.w3.org/2000/svg" viewBox="0 0 24 24" width="20" height="20" aria-hidden="true" focusable="false">< path d="M13 11.8l6.1-6.3-1-1-6.1 6.2-6.1-6.2-1 1 6.1 6.3-6.5 6.7 1 1 6.5-6.6 6.5 6.6 1-1z">

As you can see from the above chart, there is a significant difference between North America and the rest of the world. Likely, it could be related to heavy industry, considering the United States and Canada.

Now, we will be modeling the based on some variables to understand the motives behind the methane emission.

#Modeling
library(tidymodels)

#Dataset for modeling
df_mod <- 
  df_methane %>% 
  filter(region != "Aggregates",
         income != "Not classified") %>% 
  mutate(latitude = latitude %>% as.numeric(),
         longitude = longitude %>% as.numeric(),
         income = income %>% as_factor(),
         region = region %>% as_factor()) %>% 
  select(region, income, longitude,latitude, methane) 


#Splitting
set.seed(12345)
df_split <- initial_split(df_mod, 
                          prop = 0.9,
                          strata = "income")

df_train <- training(df_split)
df_test <- testing(df_split)

#10-fold cross validation for tuning
set.seed(12345)
df_folds <- 
  vfold_cv(df_train,
           strata = income,
           repeats = 5)

#Preprocessing
rec_norm <- 
  recipe(methane ~ ., data = df_train) %>% 
  step_normalize(all_numeric()) %>% 
  step_dummy(all_nominal_predictors()) 

#Model specifications

#Linear regression via brulee
library(brulee)

spec_brulee <- 
  linear_reg(penalty = tune(),
             mixture = tune()) %>% 
  set_engine("brulee")

#Multivariate adaptive regression splines (MARS)
spec_mars <- 
  mars(prod_degree = tune()) %>%
  set_engine("earth") %>% 
  set_mode("regression")

#Linear support vector machines (SVMs) via kernlab
spec_svm <- 
  svm_linear(cost = tune(),
             margin = tune()) %>%  
  set_engine("kernlab") %>% 
  set_mode("regression")

#Radial basis function support vector machines
spec_svm_rbf <- 
  svm_rbf(cost = tune(), 
          margin = tune()) %>%
  set_mode("regression") %>% 
  set_engine("kernlab")

#Workflow set
all_workflows <- 
  workflow_set(
    preproc = list(normalization = rec_norm),
    models = list(
      SVM_RBF = spec_svm_rbf,
      SVM_Linear = spec_svm,
      MARS = spec_mars,
      Brulee = spec_brulee
    )
  )

#Tuning and evaluating models
grid_ctrl <-
  control_grid(
    save_pred = TRUE,
    parallel_over = "everything",
    save_workflow = TRUE
  )

grid_results <-
  all_workflows %>%
  workflow_map(
    seed = 98765,
    resamples = df_folds,
    grid = 10,
    control = grid_ctrl
  )

#Model ranking
grid_results %>% 
  rank_results(select_best = TRUE,
               rank_metric = "rsq") %>% 
  filter(.metric == "rsq") %>% 
  select(wflow_id, rsq = mean) %>% 
  mutate(wflow_id = str_remove(wflow_id,"normalization_"))

# A tibble: 4 × 2
#  wflow_id     rsq
#  <chr>      <dbl>
#1 SVM_RBF    0.668
#2 MARS       0.625
#3 Brulee     0.176
#4 SVM_Linear 0.130

#Extracting the best model
best_results <- 
  grid_results %>% 
  extract_workflow_set_result("normalization_SVM_RBF") %>% 
  select_best(metric = "rsq")

#Finalizing the model
model_finalized <- 
  grid_results %>% 
  extract_workflow("normalization_SVM_RBF") %>% 
  finalize_workflow(best_results) %>% 
  last_fit(df_split)

#Accuracy of the finalized model
collect_metrics(model_finalized)

# A tibble: 2 × 4
#  .metric .estimator .estimate .config             
#  <chr>   <chr>          <dbl> <chr>               
#1 rmse    standard       0.534 Preprocessor1_Model1
#2 rsq     standard       0.680 Preprocessor1_Model1

As we have chosen our model we can plot the variable importance ranking.

#Variable importance
library(DALEXtra)

#Processed data frame for variable importance calculation
imp_data <- 
  rec_norm %>% 
  prep() %>% 
  bake(new_data = NULL) 

#Explainer object
explainer_svm_rbf <- 
  explain_tidymodels(
    model_finalized %>% extract_fit_parsnip(), 
    data = imp_data %>% select(-methane), 
    y = imp_data$methane,
    label = "",
    verbose = FALSE
  )


#Calculating permutation-based variable importance 
set.seed(1983)
vip_svm_rbf <- model_parts(explainer_svm_rbf, 
                           loss_function = loss_root_mean_square,
                           type = "difference",
                           B = 100,#the number of permutations
                           label = "")

#Variable importance plot
vip_svm_rbf %>% 
  mutate(variable = str_remove_all(variable, "region_|income_")) %>% 
  #removes (...) and replacing with 'and' instead
  mutate(variable = str_replace_all(variable, "\\.{3}"," and ")) %>% 
  #removes (.) and replacing with blank space
  mutate(variable = str_replace_all(variable, "\\.", " ")) %>% 
  plot() +
  labs(color = "",
       x = "",
       y = "",
       subtitle = "Higher indicates more important",
       title = "Feature Importance over 100 Permutations") +
  theme_minimal(base_family = "Bricolage Grotesque",
                base_size = 16) +
  theme(legend.position = "none",
        plot.title = element_text(hjust = 0.5, 
                                  size = 14,
                                  face = "bold"),
        plot.subtitle = element_text(hjust = 0.5, size = 12),
        panel.grid.minor.x = element_blank(),
        panel.grid.major.y = element_blank(),
        plot.background = element_rect(fill = "#F0F0F0"))

< svg xmlns="http://www.w3.org/2000/svg" width="12" height="12" fill="none" viewBox="0 0 12 12"> < path fill="#fff" d="M2 0a2 2 0 0 0-2 2v2h1.5V2a.5.5 0 0 1 .5-.5h2V0H2Zm2 10.5H2a.5.5 0 0 1-.5-.5V8H0v2a2 2 0 0 0 2 2h2v-1.5ZM8 12v-1.5h2a.5.5 0 0 0 .5-.5V8H12v2a2 2 0 0 1-2 2H8Zm2-12a2 2 0 0 1 2 2v2h-1.5V2a.5.5 0 0 0-.5-.5H8V0h2Z" />

< svg xmlns="http://www.w3.org/2000/svg" viewBox="0 0 24 24" width="20" height="20" aria-hidden="true" focusable="false">< path d="M13 11.8l6.1-6.3-1-1-6.1 6.2-6.1-6.2-1 1 6.1 6.3-6.5 6.7 1 1 6.5-6.6 6.5 6.6 1-1z">

As could be extracted from the above chart, we can say that the geospatial variables are more important than the economic levels in affecting methane gas emission.