drivers_per_season
Diego
2026-07-16
Lets create a drivers data per season
Packages
# Load packages
library(readr)
library(dplyr)## Warning: package 'dplyr' was built under R version 4.3.3##
## Attaching package: 'dplyr'## The following objects are masked from 'package:stats':
##
## filter, lag## The following objects are masked from 'package:base':
##
## intersect, setdiff, setequal, unionlibrary(ggplot2)
library(lubridate)##
## Attaching package: 'lubridate'## The following objects are masked from 'package:base':
##
## date, intersect, setdiff, unionlibrary(DataExplorer)
library(tidyverse)## ── Attaching core tidyverse packages ──────────────────────── tidyverse 2.0.0 ──
## ✔ forcats 1.0.0 ✔ tibble 3.2.1
## ✔ purrr 1.0.2 ✔ tidyr 1.3.1
## ✔ stringr 1.5.1## ── Conflicts ────────────────────────────────────────── tidyverse_conflicts() ──
## ✖ dplyr::filter() masks stats::filter()
## ✖ dplyr::lag() masks stats::lag()
## ℹ Use the conflicted package (<http://conflicted.r-lib.org/>) to force all conflicts to become errorslibrary(sf)## Warning: package 'sf' was built under R version 4.3.3## Linking to GEOS 3.11.2, GDAL 3.8.2, PROJ 9.3.1; sf_use_s2() is TRUElibrary(rnaturalearth)## Warning: package 'rnaturalearth' was built under R version 4.3.3library(lme4) ## Loading required package: Matrix
##
## Attaching package: 'Matrix'
##
## The following objects are masked from 'package:tidyr':
##
## expand, pack, unpacklibrary(emmeans) ## Welcome to emmeans.
## Caution: You lose important information if you filter this package's results.
## See '? untidy'library(tidyr)
library(broom)
library(readxl)drivers_data_mz
drivers_data_mz <- read_csv("data/processed/drivers_data.csv")## Rows: 473 Columns: 38
## ── Column specification ────────────────────────────────────────────────────────
## Delimiter: ","
## chr (14): IQR_block, IQR_TF_tubura_client, IQR_plot_fert_quality, IQR_soil_t...
## dbl (24): ICR_N_perc_23A, ICR_Org_C_avg, ICR_K_perc_23A, ICR_Avail_P_avg, IC...
##
## ℹ Use `spec()` to retrieve the full column specification for this data.
## ℹ Specify the column types or set `show_col_types = FALSE` to quiet this message.CONV_CA_Ratio <- read_csv("data/processed/CONV_CA_Ratio.csv")## Rows: 2818 Columns: 150
## ── Column specification ────────────────────────────────────────────────────────
## Delimiter: ","
## chr (49): IQR_Season, Seas_block_treat, IQR_SeasAB, Sample code, IQF_agzone...
## dbl (101): IQR_TF_tubura_client, ICR_age_hh_head, ICR_HH_size, ICR_arable_la...
##
## ℹ Use `spec()` to retrieve the full column specification for this data.
## ℹ Specify the column types or set `show_col_types = FALSE` to quiet this message.r1 <- read_csv("data/processed/r1.csv")## Rows: 2818 Columns: 33
## ── Column specification ────────────────────────────────────────────────────────
## Delimiter: ","
## chr (5): IQR_Season, IQF_District, IQR_block, cell, crop
## dbl (27): acc_rain_DBP10, acc_rain_DAP10, acc_rain_DAP20, acc_rain_DAP30, a...
## date (1): ICF_planting_date
##
## ℹ Use `spec()` to retrieve the full column specification for this data.
## ℹ Specify the column types or set `show_col_types = FALSE` to quiet this message.# vector of variables to remove
vars_to_remove <- c(
"ICR_rainfall_avg", "Comp_amount_corr_quali", "ICF_planting_date",
"avg_Y_ratio", "CA_typology", "CA_typology_09",
"acc_rain_p60_maize", "acc_rain_p60_beans",
"avg_ICR_mulch_perc_planti",
"rain_DAP60_beans", "rain_DAP120_maize",
"daily_rain_beans", "daily_rain_maize", "daily_season_rain",
"rain_DAP60_beans_chirps", "rain_DAP120_maize_chirps"
)
# remove columns
drivers_data_mz <- drivers_data_mz[ , !(names(drivers_data_mz) %in% vars_to_remove)]Now we add the season specifics: yield, compost, rain, planting date and mulch at planting
###################We add yield####################
library(dplyr)
# Step 1: compute averages per block (maize only)
avg_y_ratio_mz <- CONV_CA_Ratio %>%
filter(crop == "Maize") %>% # ⚠️ check capitalization ("Maize" vs "maize")
group_by(IQR_block) %>%
summarise(
avg_y_ratio_mz = mean(Y_ratio, na.rm = TRUE),
.groups = "drop"
)
# Step 2: join to drivers_data_mz
drivers_data_mz <- drivers_data_mz %>%
left_join(avg_y_ratio_mz, by = "IQR_block")
################## Mulch at planting###########################
library(dplyr)
# Step 1: compute average mulch % per block (maize only)
avg_mulch_mz <- CONV_CA_Ratio %>%
filter(crop == "Maize") %>% # adjust if needed ("maize")
group_by(IQR_block) %>%
summarise(
Mulch_planting_mz = mean(ICR_mulch_perc_planti, na.rm = TRUE),
.groups = "drop"
)
# Step 2: join to drivers_data_mz
drivers_data_mz <- drivers_data_mz %>%
left_join(avg_mulch_mz, by = "IQR_block")
################# COmpost ##################
#
# library(dplyr)
#
# CONV_CA_Ratio <- CONV_CA_Ratio %>%
# mutate(
# compost_quality_score = case_when(
# IQR_compost_quality == "poor" ~ 1,
# IQR_compost_quality == "average decomposed" ~ 3,
# IQR_compost_quality == "well_decomposed" ~ 5,
# TRUE ~ 2.5 # 👈 default when NA or unknown
# )
# )
# CONV_CA_Ratio <- CONV_CA_Ratio %>%
# mutate(
# Comp_amount_corr_quali =
# ICR_compost_tn_ha * ((compost_quality_score / 5) * 1.2)
# )
comp_corr_mz <- CONV_CA_Ratio %>%
filter(crop == "Maize") %>% # adjust if needed
group_by(IQR_block) %>%
summarise(
Comp_amount_corr_quali_mz = mean(Comp_amount_corr_quali, na.rm = TRUE),
.groups = "drop"
)
drivers_data_mz <- drivers_data_mz %>%
left_join(comp_corr_mz, by = "IQR_block")
##################Planting date################
# library(dplyr)
# library(lubridate)
# library(stringr)
#
# ###########################################################
# ### 1. Parse both date formats + detect season type
# ###########################################################
#
# CONV_CA_Ratio <- CONV_CA_Ratio %>%
# mutate(
# planting_dmy = suppressWarnings(dmy(ICF_planting_date)),
# planting_mdy = suppressWarnings(mdy(ICF_planting_date)),
# season_type = case_when(
# str_detect(IQR_Season, "A$") ~ "A",
# str_detect(IQR_Season, "B$") ~ "B",
# TRUE ~ NA_character_
# )
# )
#
# ###########################################################
# ### 2. Expected planting months
# ###########################################################
#
# expected_months <- list(
# A = 8:10, # Aug–Oct
# B = 1:3 # Jan–Mar
# )
#
# ###########################################################
# ### 3. Decide correct format per season × district
# ###########################################################
#
# format_decision <- CONV_CA_Ratio %>%
# filter(!is.na(IQR_Season), !is.na(IQF_District)) %>%
# group_by(IQR_Season, IQF_District, season_type) %>%
# summarise(
# dmy_good = sum(month(planting_dmy) %in% expected_months[[first(season_type)]], na.rm = TRUE),
# mdy_good = sum(month(planting_mdy) %in% expected_months[[first(season_type)]], na.rm = TRUE),
# best_format = case_when(
# dmy_good > mdy_good ~ "dmy",
# mdy_good > dmy_good ~ "mdy",
# TRUE ~ "unknown"
# ),
# .groups = "drop"
# )
#
# ###########################################################
# ### 4. Join and clean planting date
# ###########################################################
#
# CONV_CA_Ratio <- CONV_CA_Ratio %>%
# select(-any_of("best_format")) %>%
# left_join(
# format_decision %>%
# select(IQR_Season, IQF_District, best_format),
# by = c("IQR_Season", "IQF_District")
# ) %>%
# mutate(
# ICF_planting_date_clean = case_when(
# best_format == "dmy" ~ planting_dmy,
# best_format == "mdy" ~ planting_mdy,
# TRUE ~ NA_Date_
# )
# )
#
# ###########################################################
# ### 5. Compute planting DOY
# ###########################################################
#
# CONV_CA_Ratio <- CONV_CA_Ratio %>%
# mutate(planting_doy = yday(ICF_planting_date_clean))
#
# ###########################################################
# ### 6. Earliest planting per district-season
# ###########################################################
#
# earliest_doy_by_district_season <- CONV_CA_Ratio %>%
# filter(!is.na(planting_doy)) %>%
# group_by(IQR_Season, IQF_District) %>%
# summarise(first_doy = min(planting_doy), .groups = "drop")
#
# ###########################################################
# ### 7. Compute planting delay
# ###########################################################
#
# CONV_CA_Ratio <- CONV_CA_Ratio %>%
# select(-starts_with("first_doy")) %>%
# left_join(
# earliest_doy_by_district_season,
# by = c("IQR_Season", "IQF_District")
# ) %>%
# mutate(
# planting_delay_doy = planting_doy - first_doy
# )
# Step 1: average planting delay per block (maize only)
planting_delay_mz <- CONV_CA_Ratio %>%
filter(crop == "Maize") %>% # adjust if needed
group_by(IQR_block) %>%
summarise(
Planting_date_rel_mz = mean(planting_delay_doy, na.rm = TRUE),
.groups = "drop"
)
# Step 2: join to drivers_data_mz
drivers_data_mz <- drivers_data_mz %>%
left_join(planting_delay_mz, by = "IQR_block")
############################Rain################################
# accumulated rain 0-30, 0-60 amd and 60-12 0-120 for maize
library(dplyr)
avg_rain_mz <- r1 %>%
filter(crop == "Maize") %>% # ⚠️ adjust if "maize"
group_by(IQR_block) %>%
summarise(
acc_rain_DAP30_mz = mean(acc_rain_DAP30, na.rm = TRUE),
acc_rain_DAP60_mz = mean(acc_rain_DAP60, na.rm = TRUE),
acc_rain_DAP120_mz = mean(acc_rain_DAP120, na.rm = TRUE),
acc_rain_DAP60_120_mz = mean(acc_rain_DAP60_120, na.rm = TRUE),
.groups = "drop"
)
drivers_data_mz <- drivers_data_mz %>%
left_join(avg_rain_mz, by = "IQR_block")| ## Now for maize but only 24 adn 25 #```{r} ###################We add yield#################### library(dplyr) |
| # Step 1: compute averages per block (maize only) avg_y_ratio_mz <- CONV_CA_Ratio %>% filter(IQR_Season %in% c(“24A”, “25A”)) %>% # ⚠️ check capitalization (“Maize” vs “maize”) group_by(IQR_block) %>% summarise( avg_y_ratio_mz = mean(Y_ratio, na.rm = TRUE), .groups = “drop” ) |
| # Step 2: join to drivers_data_mz drivers_data_mz <- drivers_data_mz %>% left_join(avg_y_ratio_mz, by = “IQR_block”) |
| ################## Mulch at planting########################### library(dplyr) |
| # Step 1: compute average mulch % per block (maize only) avg_mulch_mz <- CONV_CA_Ratio %>% filter(IQR_Season %in% c(“24A”, “25A”)) %>% # adjust if needed (“maize”) group_by(IQR_block) %>% summarise( Mulch_planting_mz = mean(ICR_mulch_perc_planti, na.rm = TRUE), .groups = “drop” ) |
| # Step 2: join to drivers_data_mz drivers_data_mz <- drivers_data_mz %>% left_join(avg_mulch_mz, by = “IQR_block”) |
| ################# COmpost ################## |
| comp_corr_mz <- CONV_CA_Ratio %>% filter(IQR_Season %in% c(“24A”, “25A”)) %>% # adjust if needed group_by(IQR_block) %>% summarise( Comp_amount_corr_quali_mz = mean(Comp_amount_corr_quali, na.rm = TRUE), .groups = “drop” ) |
| drivers_data_mz <- drivers_data_mz %>% left_join(comp_corr_mz, by = “IQR_block”) |
| # Step 1: average planting delay per block (maize only) planting_delay_mz <- CONV_CA_Ratio %>% filter(IQR_Season %in% c(“24A”, “25A”)) %>% # adjust if needed group_by(IQR_block) %>% summarise( Planting_date_rel_mz = mean(planting_delay_doy, na.rm = TRUE), .groups = “drop” ) |
| # Step 2: join to drivers_data_mz drivers_data_mz <- drivers_data_mz %>% left_join(planting_delay_mz, by = “IQR_block”) |
| ############################Rain################################ |
| # accumulated rain 0-30, 0-60 amd and 60-12 0-120 for maize |
| library(dplyr) |
| avg_rain_mz <- r1 %>% filter(IQR_Season %in% c(“24A”, “25A”))%>% # ⚠️ adjust if “maize” group_by(IQR_block) %>% summarise( acc_rain_DAP30_mz = mean(acc_rain_DAP30, na.rm = TRUE), acc_rain_DAP60_mz = mean(acc_rain_DAP60, na.rm = TRUE), acc_rain_DAP120_mz = mean(acc_rain_DAP120, na.rm = TRUE), acc_rain_DAP60_120_mz = mean(acc_rain_DAP60_120, na.rm = TRUE), .groups = “drop” ) drivers_data_mz <- drivers_data_mz %>% left_join(avg_rain_mz, by = “IQR_block”) |
| #``` |
| ## drivers_data_bn |
r drivers_data_bn <- read_csv("data/processed/drivers_data.csv") |
## Rows: 473 Columns: 38 ## ── Column specification ──────────────────────────────────────────────────────── ## Delimiter: "," ## chr (14): IQR_block, IQR_TF_tubura_client, IQR_plot_fert_quality, IQR_soil_t... ## dbl (24): ICR_N_perc_23A, ICR_Org_C_avg, ICR_K_perc_23A, ICR_Avail_P_avg, IC... ## ## ℹ Use `spec()` to retrieve the full column specification for this data. ## ℹ Specify the column types or set `show_col_types = FALSE` to quiet this message. |
| ```r # vector of variables to remove vars_to_remove <- c( “ICR_rainfall_avg”, “Comp_amount_corr_quali”, “ICF_planting_date”, “avg_Y_ratio”, “CA_typology”, “CA_typology_09”, “acc_rain_p60_maize”, “acc_rain_p60_beans”, “avg_ICR_mulch_perc_planti”, “rain_DAP60_beans”, “rain_DAP120_maize”, “daily_rain_beans”, “daily_rain_maize”, “daily_season_rain”, “rain_DAP60_beans_chirps”, “rain_DAP120_maize_chirps” ) |
| # remove columns drivers_data_bn <- drivers_data_bn[ , !(names(drivers_data_bn) %in% vars_to_remove)] ``` Now we add the season specifics: yield, compost, rain, planting date and mulch at planting |
| ```r ###################We add yield#################### library(dplyr) |
| # Step 1: compute averages per block (maize only) avg_y_ratio_bn <- CONV_CA_Ratio %>% filter(crop == “Beans”) %>% # ⚠️ check capitalization (“Maize” vs “maize”) group_by(IQR_block) %>% summarise( avg_y_ratio_bn = mean(Y_ratio, na.rm = TRUE), .groups = “drop” ) |
| # Step 2: join to drivers_data_bn drivers_data_bn <- drivers_data_bn %>% left_join(avg_y_ratio_bn, by = “IQR_block”) |
| ################## Mulch at planting########################### library(dplyr) |
| # Step 1: compute average mulch % per block (maize only) avg_mulch_bn <- CONV_CA_Ratio %>% filter(crop == “Beans”) %>% # adjust if needed (“maize”) group_by(IQR_block) %>% summarise( Mulch_planting_bn = mean(ICR_mulch_perc_planti, na.rm = TRUE), .groups = “drop” ) |
| # Step 2: join to drivers_data_mz drivers_data_bn <- drivers_data_bn %>% left_join(avg_mulch_bn, by = “IQR_block”) |
| ################# COmpost ################## |
| comp_corr_bn <- CONV_CA_Ratio %>% filter(crop == “Beans”) %>% # adjust if needed group_by(IQR_block) %>% summarise( Comp_amount_corr_quali_bn = mean(Comp_amount_corr_quali, na.rm = TRUE), .groups = “drop” ) |
| drivers_data_bn <- drivers_data_bn %>% left_join(comp_corr_bn, by = “IQR_block”) |
| ##################Planting date################ |
| # Step 1: average planting delay per block (maize only) planting_delay_bn <- CONV_CA_Ratio %>% filter(crop == “Beans”) %>% # adjust if needed group_by(IQR_block) %>% summarise( Planting_date_rel_bn = mean(planting_delay_doy, na.rm = TRUE), .groups = “drop” ) |
| # Step 2: join to drivers_data_mz drivers_data_bn <- drivers_data_bn %>% left_join(planting_delay_bn, by = “IQR_block”) |
| ############################Rain################################ |
| # accumulated rain 0-30, 0-60 amd and 60-12 0-120 for maize |
| library(dplyr) |
| avg_rain_bn <- r1 %>% filter(crop == “Beans”) %>% # ⚠️ adjust if “maize” group_by(IQR_block) %>% summarise( acc_rain_DAP30_bn = mean(acc_rain_DAP30, na.rm = TRUE), acc_rain_DAP60_bn = mean(acc_rain_DAP60, na.rm = TRUE), acc_rain_DAP120_bn = mean(acc_rain_DAP120, na.rm = TRUE), acc_rain_DAP60_120_bn = mean(acc_rain_DAP60_120, na.rm = TRUE), .groups = “drop” ) drivers_data_bn <- drivers_data_bn %>% left_join(avg_rain_bn, by = “IQR_block”) ``` |
| ## Now for beans but only 24 adn 25 #```{r} ###################We add yield#################### library(dplyr) |
| # Step 1: compute averages per block (maize only) avg_y_ratio_bn <- CONV_CA_Ratio %>% filter(IQR_Season %in% c(“24A”, “25A”)) %>% # ⚠️ check capitalization (“Maize” vs “maize”) group_by(IQR_block) %>% summarise( avg_y_ratio_bn = mean(Y_ratio, na.rm = TRUE), .groups = “drop” ) |
| # Step 2: join to drivers_data_mz drivers_data_bn <- drivers_data_bn %>% left_join(avg_y_ratio_bn, by = “IQR_block”) |
| ################## Mulch at planting########################### library(dplyr) |
| # Step 1: compute average mulch % per block (maize only) avg_mulch_bn <- CONV_CA_Ratio %>% filter(IQR_Season %in% c(“24A”, “25A”)) %>% # adjust if needed (“maize”) group_by(IQR_block) %>% summarise( Mulch_planting_bn = mean(ICR_mulch_perc_planti, na.rm = TRUE), .groups = “drop” ) |
| # Step 2: join to drivers_data_mz drivers_data_bn <- drivers_data_bn %>% left_join(avg_mulch_bn, by = “IQR_block”) |
| ################# COmpost ################## |
| comp_corr_bn <- CONV_CA_Ratio %>% filter(IQR_Season %in% c(“24A”, “25A”)) %>% # adjust if needed group_by(IQR_block) %>% summarise( Comp_amount_corr_quali_bn = mean(Comp_amount_corr_quali, na.rm = TRUE), .groups = “drop” ) |
| drivers_data_bn <- drivers_data_bn %>% left_join(comp_corr_bn, by = “IQR_block”) |
| # Step 1: average planting delay per block (maize only) planting_delay_bn <- CONV_CA_Ratio %>% filter(IQR_Season %in% c(“24A”, “25A”)) %>% # adjust if needed group_by(IQR_block) %>% summarise( Planting_date_rel_bn = mean(planting_delay_doy, na.rm = TRUE), .groups = “drop” ) |
| # Step 2: join to drivers_data_mz drivers_data_bn <- drivers_data_bn %>% left_join(planting_delay_bn, by = “IQR_block”) |
| ############################Rain################################ |
| # accumulated rain 0-30, 0-60 amd and 60-12 0-120 for maize |
| library(dplyr) |
| avg_rain_bn <- r1 %>% filter(IQR_Season %in% c(“24A”, “25A”))%>% # ⚠️ adjust if “maize” group_by(IQR_block) %>% summarise( acc_rain_DAP30_bn = mean(acc_rain_DAP30, na.rm = TRUE), acc_rain_DAP60_bn = mean(acc_rain_DAP60, na.rm = TRUE), acc_rain_DAP120_bn = mean(acc_rain_DAP120, na.rm = TRUE), acc_rain_DAP60_120_bn = mean(acc_rain_DAP60_120, na.rm = TRUE), .groups = “drop” ) drivers_data_bn <- drivers_data_bn %>% left_join(avg_rain_bn, by = “IQR_block”) |
| #``` |
| # Random Forest for Maize |
| # Lets firs check for outlier for y ration for maize |
| ```r library(dplyr) |
| drivers_data_mz <- drivers_data_mz %>% group_by(IQF_agzone) %>% mutate( Q1 = quantile(avg_y_ratio_mz, 0.25, na.rm = TRUE), Q3 = quantile(avg_y_ratio_mz, 0.75, na.rm = TRUE), IQR_val = Q3 - Q1, lower_bound = Q1 - 1.5 * IQR_val, upper_bound = Q3 + 1.5 * IQR_val, outlier_y_ratio = avg_y_ratio_mz < lower_bound | avg_y_ratio_mz > upper_bound ) %>% ungroup() |
| drivers_data_mz <- drivers_data_mz %>% filter( !is.na(avg_y_ratio_mz), # remove missing values !outlier_y_ratio # remove outliers ) ``` |
| ## First lets check if the is NA and what to do with them |
r colSums(is.na(drivers_data_mz)) |
## IQR_block ICR_N_perc_23A ICR_Org_C_avg ## 0 11 0 ## ICR_K_perc_23A ICR_Avail_P_avg ICR_arable_land_avg ## 11 0 0 ## ICF_Alt_avg IQR_TF_tubura_client IQR_plot_fert_quality ## 0 0 2 ## IQR_soil_texture IQR_soil_color IQF_position_slope ## 0 0 0 ## Weed_species_A_combined IQF_environment Slope ## 0 0 0 ## K_factor P_factor slope_cat ## 11 0 0 ## Couch_Consit Weed_Overall_avg Whitch_weed ## 0 0 0 ## IQF_agzone avg_y_ratio_mz Mulch_planting_mz ## 0 0 0 ## Comp_amount_corr_quali_mz Planting_date_rel_mz acc_rain_DAP30_mz ## 0 0 0 ## acc_rain_DAP60_mz acc_rain_DAP120_mz acc_rain_DAP60_120_mz ## 0 0 0 ## Q1 Q3 IQR_val ## 0 0 0 ## lower_bound upper_bound outlier_y_ratio ## 0 0 0
# we complete them with average of all the other int he same AEZ |
| ```r library(dplyr) |
| drivers_data_mz_c <- drivers_data_mz %>% group_by(IQF_agzone) %>% mutate( across( where(is.numeric), ~ ifelse(is.na(.), mean(., na.rm = TRUE), .) ) ) %>% ungroup() ``` |
| ```r library(readr) library(dplyr) library(ggplot2) library(lubridate) library(DataExplorer) library(tidyverse) library(sf) library(rnaturalearth) library(lme4) library(emmeans) library(tidyr) library(broom) library(readxl) |
| RF_data_mz <- drivers_data_mz_c |
| # Convert binary variables (0/1) to factors RF_data_mz\(Couch_Consit <- factor(RF_data_mz\)Couch_Consit, levels = c(0,1), labels = c(“Absent”,“Present”)) |
| RF_data_mz\(Whitch_weed <- factor(RF_data_mz\)Whitch_weed, levels = c(0,1), labels = c(“Absent”,“Present”)) |
| # Convert categorical variables to factors RF_data_mz\(IQR_soil_texture <- as.factor(RF_data_mz\)IQR_soil_texture) RF_data_mz\(IQR_soil_color <- as.factor(RF_data_mz\)IQR_soil_color) RF_data_mz\(IQF_agzone <- as.factor(RF_data_mz\)IQF_agzone) RF_data_mz\(IQF_environment <- as.factor(RF_data_mz\)IQF_environment) RF_data_mz\(IQF_position_slope <- as.factor(RF_data_mz\)IQF_position_slope) |
| # Optional (often categorical as well) RF_data_mz\(IQR_plot_fert_quality <- as.factor(RF_data_mz\)IQR_plot_fert_quality) ``` |
| # Select variables and prepare dataset |
| lets first check which of the rains have more correlation with Yratio |
| ```r vars <- c( “acc_rain_DAP120_mz”, “acc_rain_DAP30_mz”, “acc_rain_DAP60_mz”, “acc_rain_DAP60_120_mz” ) |
| results <- lapply(vars, function(v) { formula <- as.formula(paste0(“avg_y_ratio_mz ~”, v, ” + I(“, v,”^2)“)) model <- lm(formula, data = RF_data_mz) |
| summary_model <- summary(model) |
| data.frame( variable = v, R2 = summary_model\(r.squared, p_linear = summary_model\)coefficients[2, 4], p_quadratic = summary_model$coefficients[3, 4] ) }) |
| quad_results <- do.call(rbind, results) quad_results ``` |
## variable R2 p_linear p_quadratic ## 1 acc_rain_DAP120_mz 0.02146602 0.003041694 0.001760938 ## 2 acc_rain_DAP30_mz 0.01200727 0.024200070 0.019531963 ## 3 acc_rain_DAP60_mz 0.01699926 0.009169896 0.005484931 ## 4 acc_rain_DAP60_120_mz 0.01427346 0.015211793 0.010864291 |
r library(ggplot2) library(dplyr) library(gridExtra) |
## Warning: package 'gridExtra' was built under R version 4.3.3 |
## ## Attaching package: 'gridExtra' |
## The following object is masked from 'package:dplyr': ## ## combine |
| ```r vars <- c( “acc_rain_DAP120_mz”, “acc_rain_DAP30_mz”, “acc_rain_DAP60_mz”, “acc_rain_DAP60_120_mz” ) |
| plots <- lapply(vars, function(v) { |
| # Build formula formula <- as.formula(paste0(“avg_y_ratio_mz ~”, v, ” + I(“, v,”^2)“)) |
| # Fit model model <- lm(formula, data = RF_data_mz) summary_model <- summary(model) |
| # Extract stats r2 <- round(summary_model\(r.squared, 3) p_lin <- signif(summary_model\)coefficients[2, 4], 2) p_quad <- signif(summary_model$coefficients[3, 4], 2) |
| # Create plot ggplot(RF_data_mz, aes_string(x = v, y = “avg_y_ratio_mz”)) + geom_point(alpha = 0.6) + stat_smooth(method = “lm”, formula = y ~ x + I(x^2), color = “blue”) + theme_minimal() + labs( title = v, subtitle = paste0(“R² =”, r2, ” | p(linear) = “, p_lin,” | p(quadratic) = “, p_quad), x =”Rain”, y = “avg_y_ratio_mz” ) }) ``` |
## Warning: `aes_string()` was deprecated in ggplot2 3.0.0. ## ℹ Please use tidy evaluation idioms with `aes()`. ## ℹ See also `vignette("ggplot2-in-packages")` for more information. ## This warning is displayed once every 8 hours. ## Call `lifecycle::last_lifecycle_warnings()` to see where this warning was ## generated. |
r # Arrange all plots together grid.arrange(grobs = plots, ncol = 2) |
 |
r library(tidyverse) library(caret) |
## Warning: package 'caret' was built under R version 4.3.3 |
## Loading required package: lattice |
## ## Attaching package: 'caret' |
## The following object is masked from 'package:purrr': ## ## lift |
r library(randomForest) |
## Warning: package 'randomForest' was built under R version 4.3.3 |
## randomForest 4.7-1.2 |
## Type rfNews() to see new features/changes/bug fixes. |
## ## Attaching package: 'randomForest' |
## The following object is masked from 'package:gridExtra': ## ## combine |
## The following object is masked from 'package:ggplot2': ## ## margin |
## The following object is masked from 'package:dplyr': ## ## combine |
r library(usdm) |
## Warning: package 'usdm' was built under R version 4.3.3 |
## Loading required package: terra |
## Warning: package 'terra' was built under R version 4.3.3 |
## terra 1.8.29 |
## ## Attaching package: 'terra' |
## The following object is masked from 'package:tidyr': ## ## extract |
r # Select response + explanatory variables data_model <- RF_data_mz %>% dplyr::select( avg_y_ratio_mz, ICR_N_perc_23A, #ICR_K_perc_23A, K_factor, Couch_Consit, Whitch_weed, ICF_Alt_avg, #ICR_Avail_P_avg, acc_rain_DAP120_mz, acc_rain_DAP30_mz, #acc_rain_DAP60_mz, #acc_rain_DAP60_120_mz, P_factor, #ICR_arable_land_avg, Comp_amount_corr_quali_mz, Planting_date_rel_mz, Slope, #IQR_TF_tubura_client, IQR_plot_fert_quality, IQR_soil_texture, ICR_Org_C_avg, IQR_soil_color, IQF_position_slope, Weed_species_A_combined, Weed_Overall_avg, Mulch_planting_mz, #IQF_agzone, #IQF_environment ) %>% na.omit() |
| # Train/ Test split |
| ```r set.seed(42) |
| parts <- createDataPartition(data_model$avg_y_ratio_mz, p = 0.7, list = FALSE) |
| data_train <- data_model[parts, ] data_test <- data_model[-parts, ] |
| x_train <- data_train[, -1] y_train <- data_train$avg_y_ratio_mz ``` |
| # Variable selection (RFE) |
| ```r set.seed(627) |
| control_RFE <- rfeControl( functions = rfFuncs, method = “repeatedcv”, repeats = 1, number = 5 ) |
| start <- Sys.time() |
| result_RFE <- rfe( x = x_train, y = y_train, sizes = 1:ncol(x_train), rfeControl = control_RFE ) |
| end <- Sys.time() print(end - start) ``` |
## Time difference of 29.09347 secs |
| # Selected: |
r predictors(result_RFE) |
## [1] "Weed_Overall_avg" "ICF_Alt_avg" ## [3] "Couch_Consit" "Comp_amount_corr_quali_mz" ## [5] "acc_rain_DAP120_mz" |
| # non selected |
r names(x_train)[!names(x_train) %in% result_RFE$optVariables] |
## [1] "ICR_N_perc_23A" "K_factor" ## [3] "Whitch_weed" "acc_rain_DAP30_mz" ## [5] "P_factor" "Planting_date_rel_mz" ## [7] "Slope" "IQR_plot_fert_quality" ## [9] "IQR_soil_texture" "ICR_Org_C_avg" ## [11] "IQR_soil_color" "IQF_position_slope" ## [13] "Weed_species_A_combined" "Mulch_planting_mz" |
| # Keep only selected variables |
r dataRFE <- subset(data_train, select = result_RFE$optVariables) |
| I will force it to keep “Slope”, “IQR_soil_texture”, “acc_rain_p60_maize”, “ICR_Avail_P_avg”, “Comp_amount_corr_quali”, “IQF_position_slope”, “IQR_soil_texture”, “IQR_soil_color” |
| ```r forced_vars <- c(“Mulch_planting_mz”,“Planting_date_rel_mz”,“acc_rain_DAP30_mz”, “Slope”,“IQR_soil_color”,“IQF_position_slope”, “K_factor”,“ICR_Org_C_avg” , “P_factor”, “ICR_N_perc_23A” ,“IQR_soil_texture” ) |
| final_vars <- unique(c(result_RFE$optVariables, forced_vars)) |
| dataRFE <- data_train[, final_vars] ``` |
| # Test colinearity |
| ## First we separate numeric from factoria as VIF only handles numeric |
| ```r library(dplyr) |
| data_numeric <- dataRFE %>% dplyr::select(where(is.numeric)) data_factor <- dataRFE %>% dplyr::select(where(is.factor)) |
| names(data_numeric) ``` |
## [1] "Weed_Overall_avg" "ICF_Alt_avg" ## [3] "Comp_amount_corr_quali_mz" "acc_rain_DAP120_mz" ## [5] "Mulch_planting_mz" "Planting_date_rel_mz" ## [7] "acc_rain_DAP30_mz" "Slope" ## [9] "ICR_Org_C_avg" "ICR_N_perc_23A" |
r names(data_factor) |
## [1] "Couch_Consit" "IQR_soil_color" "IQF_position_slope" ## [4] "IQR_soil_texture" |
| ## Colineary test: VIF Common thresholds: |
| VIF < 5 → safe VIF < 10 → acceptable VIF > 10 → problematic collinearity |
| So we keep all the numeric |
| ```r library(usdm) data_numeric_df <- as.data.frame(data_numeric) |
| vif_step <- vifstep(data_numeric_df, th = 10) |
| vif_step ``` |
## No variable from the 10 input variables has collinearity problem. ## ## The linear correlation coefficients ranges between: ## min correlation ( ICR_N_perc_23A ~ Planting_date_rel_mz ): -0.0268504 ## max correlation ( acc_rain_DAP30_mz ~ acc_rain_DAP120_mz ): 0.8361558 ## ## ---------- VIFs of the remained variables -------- ## Variables VIF ## 1 Weed_Overall_avg 1.272723 ## 2 ICF_Alt_avg 1.795888 ## 3 Comp_amount_corr_quali_mz 1.113908 ## 4 acc_rain_DAP120_mz 3.742046 ## 5 Mulch_planting_mz 1.321458 ## 6 Planting_date_rel_mz 1.296388 ## 7 acc_rain_DAP30_mz 4.008136 ## 8 Slope 1.524420 ## 9 ICR_Org_C_avg 1.850698 ## 10 ICR_N_perc_23A 1.790668 |
| # Final RF database |
r data_train_final <- cbind( avg_y_ratio = data_train$avg_y_ratio_mz, dataRFE ) |
| # RF model calibration |
| ```r # model calibration ——————————————————- |
| fitControl <- trainControl( method = “repeatedcv”, repeats = 1, number = 5, search = “grid” ) |
| # grid for ranger RF rfGrid <- expand.grid( mtry = c(1:length(dataRFE)), min.node.size = c(5, 10, 50, 100), splitrule = “variance” ) |
| set.seed(627) |
| start <- Sys.time() |
| RF_model_train <- train( formula(paste0(“avg_y_ratio_mz ~”, paste0(names(dataRFE), collapse = ” + “))), data = data_train, method =”ranger”, trControl = fitControl, tuneGrid = rfGrid, metric = “RMSE”, verbose = FALSE ) |
| end <- Sys.time() print(end - start) ``` |
## Time difference of 22.79162 secs |
| # Final model using best hyperparameters |
| ```r fitControl <- trainControl(method = “none”, search = “grid”) |
| RFFGrid <- expand.grid( mtry = RF_model_train\(bestTune\)mtry, splitrule = RF_model_train\(bestTune\)splitrule, min.node.size = RF_model_train\(bestTune\)min.node.size ) |
| RFFGrid ``` |
## mtry splitrule min.node.size ## 1 3 variance 10 |
| # Train final RF |
r fit_RF <- train( formula(paste0("avg_y_ratio_mz ~ ", paste0(names(dataRFE), collapse = " + "))), data = data_train, method = "ranger", trControl = fitControl, metric = "RMSE", verbose = FALSE, tuneGrid = RFFGrid ) |
| # Train set |
r library(Metrics) |
## Warning: package 'Metrics' was built under R version 4.3.3 |
## ## Attaching package: 'Metrics' |
## The following objects are masked from 'package:caret': ## ## precision, recall |
| ```r results_RF_train <- data.frame( avg_y_ratio_mz = data_train$avg_y_ratio_mz, y_pred = predict(fit_RF, newdata = data_train) ) |
| results_RF_train\(residual <- results_RF_train\)avg_y_ratio_mz - results_RF_train$y_pred |
| RMSE_RF_train <- Metrics::rmse(results_RF_train\(avg_y_ratio_mz, results_RF_train\)y_pred) RMSE_RF_train ``` |
## [1] 0.1044726 |
r R2_RF_train <- cor(results_RF_train$y_pred, results_RF_train$avg_y_ratio_mz)^2 R2_RF_train |
## [1] 0.8326022 |
| # For the test set: |
| ```r results_RF_test <- data.frame( avg_y_ratio_mz = data_test$avg_y_ratio_mz, y_pred = predict(fit_RF, newdata = data_test) ) |
| results_RF_test\(residual <- results_RF_test\)avg_y_ratio_mz - results_RF_test$y_pred |
| RMSE_RF_test <- Metrics::rmse(results_RF_test\(avg_y_ratio_mz, results_RF_test\)y_pred) RMSE_RF_test ``` |
## [1] 0.1469754 |
r R2_RF_test <- cor(results_RF_test$y_pred, results_RF_test$avg_y_ratio_mz)^2 R2_RF_test |
## [1] 0.3570728 |
| # Observed vs Predicted plots ## Trained data |
| ```r plot(results_RF_train\(y_pred, results_RF_train\)avg_y_ratio_mz, xlab = “Predicted Y ratio”, ylab = “Observed Y ratio”, pch = 20, main = “RF Train”) |
| abline(lm(avg_y_ratio_mz ~ y_pred, data = results_RF_train), col = “red”) abline(0,1,col=“red”, lty=2) |
| legend(“topleft”, legend = c(paste(“R2 =”, round(R2_RF_train,3)), paste(“RMSE =”, round(RMSE_RF_train,3))), bty=“n”) ``` |
 |
| ## tested |
| ```r # Counts true_good <- sum(results_RF_test\(avg_y_ratio_mz >= 0.9 & results_RF_test\)y_pred >= 0.9) false_acc <- sum(results_RF_test\(avg_y_ratio_mz < 0.9 & results_RF_test\)y_pred >= 0.9) |
| plot(results_RF_test\(y_pred, results_RF_test\)avg_y_ratio_mz, xlab = “Predicted CA:CT yield ratio”, ylab = “Observed CA:CT yield ratio”, pch = 20, main = “RF Test”) |
| abline(lm(avg_y_ratio_mz ~ y_pred, data = results_RF_test), col = “red”) abline(0,1,col=“red”, lty=2) |
| # Threshold lines abline(v = 0.9, col = “blue”, lty = 2) abline(h = 0.9, col = “blue”, lty = 2) |
| # Add quadrant labels text(x = 0.895, y = 1.1, labels = paste(“True performer =”, true_good), col = “darkgreen”, pos = 4) |
| text(x = 0.895, y = 0.6, labels = paste(“False performer =”, false_acc), col = “orange”, pos = 4) |
| legend(“topleft”, legend = c(paste(“R² =”, round(R2_RF_test,3)), paste(“RMSE =”, round(RMSE_RF_test,3))), bty=“n”) ``` |
 |
| ## Stave model objects |
r save( data_model, parts, data_train, data_test, result_RFE, RF_model_train, fit_RF, results_RF_train, RMSE_RF_train, R2_RF_train, results_RF_test, RMSE_RF_test, R2_RF_test, file = "RF_avg_Y_ratio_model.Rdata" )
## names(dataRFE) %in% names(data_train) |
r names(dataRFE) %in% names(data_train) |
## [1] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE ## [16] TRUE
# Visualizetion of variables importance |
| ```r # Visualizetion of variables importance |
| # Visualizetion of variables importance |
| library(dplyr) library(stringr) library(ggplot2) |
| fit_RF <- train( formula(paste0(“avg_y_ratio_mz ~”, paste0(names(dataRFE), collapse = ” + “))), data = data_train, method =”ranger”, trControl = fitControl, metric = “RMSE”, verbose = FALSE, tuneGrid = RFFGrid, importance = “permutation” ) |
| # Extract importance importance_RF <- varImp(fit_RF) |
| # —————————– # 🔽 ADD THIS PART HERE # —————————– imp <- importance_RF\(importance imp\)variable <- rownames(imp) |
| imp_grouped <- imp %>% mutate( base_var = case_when( str_detect(variable, “^K_factor”) ~ “K_factor”, str_detect(variable, “^P_factor”) ~ “P_factor”, # ADD THIS str_detect(variable, “^Couch_Consit”) ~ “Couch_Consit”, str_detect(variable, “^Whitch_weed”) ~ “Whitch_weed”, str_detect(variable, “^IQR_soil_color”) ~ “IQR_soil_color”, str_detect(variable, “^IQR_soil_texture”) ~ “IQR_soil_texture”, str_detect(variable, “^IQF_position_slope”) ~ “IQF_position_slope”, str_detect(variable, “^IQR_plot_fert_quality”) ~ “IQR_plot_fert_quality”, str_detect(variable, “^Weed_species_A_combined”) ~ “Weed_species_A_combined”, # ADD THIS TRUE ~ variable ) ) %>% group_by(base_var) %>% summarise( Overall = max(Overall), .groups = “drop” ) %>% arrange(desc(Overall)) |
| # —————————– # Plot aggregated importance # —————————– ggplot(imp_grouped, aes(x = reorder(base_var, Overall), y = Overall)) + geom_col() + coord_flip() + theme_minimal() ``` |
 |
#NOw RF for beans
Lets firs check for outlier for y ration for maize
library(dplyr)
drivers_data_bn <- drivers_data_bn %>%
group_by(IQF_agzone) %>%
mutate(
Q1 = quantile(avg_y_ratio_bn, 0.25, na.rm = TRUE),
Q3 = quantile(avg_y_ratio_bn, 0.75, na.rm = TRUE),
IQR_val = Q3 - Q1,
lower_bound = Q1 - 1.5 * IQR_val,
upper_bound = Q3 + 1.5 * IQR_val,
outlier_y_ratio = avg_y_ratio_bn < lower_bound | avg_y_ratio_bn > upper_bound
) %>%
ungroup()
drivers_data_bn <- drivers_data_bn %>%
filter(
!is.na(avg_y_ratio_bn), # remove missing values
!outlier_y_ratio # remove outliers
)First lets check if the is NA and what to do with them
colSums(is.na(drivers_data_bn))## IQR_block ICR_N_perc_23A ICR_Org_C_avg
## 0 13 0
## ICR_K_perc_23A ICR_Avail_P_avg ICR_arable_land_avg
## 13 0 0
## ICF_Alt_avg IQR_TF_tubura_client IQR_plot_fert_quality
## 0 0 2
## IQR_soil_texture IQR_soil_color IQF_position_slope
## 0 0 0
## Weed_species_A_combined IQF_environment Slope
## 0 0 0
## K_factor P_factor slope_cat
## 13 0 0
## Couch_Consit Weed_Overall_avg Whitch_weed
## 0 0 0
## IQF_agzone avg_y_ratio_bn Mulch_planting_bn
## 0 0 0
## Comp_amount_corr_quali_bn Planting_date_rel_bn acc_rain_DAP30_bn
## 0 2 0
## acc_rain_DAP60_bn acc_rain_DAP120_bn acc_rain_DAP60_120_bn
## 0 0 0
## Q1 Q3 IQR_val
## 0 0 0
## lower_bound upper_bound outlier_y_ratio
## 0 0 0we complete them with average of all the other int he same AEZ
library(dplyr)
drivers_data_bn_c <- drivers_data_bn %>%
group_by(IQF_agzone) %>%
mutate(
across(
where(is.numeric),
~ ifelse(is.na(.), mean(., na.rm = TRUE), .)
)
) %>%
ungroup()library(readr)
library(dplyr)
library(ggplot2)
library(lubridate)
library(DataExplorer)
library(tidyverse)
library(sf)
library(rnaturalearth)
library(lme4)
library(emmeans)
library(tidyr)
library(broom)
library(readxl)
RF_data_bn <- drivers_data_bn_c
# Convert binary variables (0/1) to factors
RF_data_bn$Couch_Consit <- factor(RF_data_bn$Couch_Consit,
levels = c(0,1),
labels = c("Absent","Present"))
RF_data_bn$Whitch_weed <- factor(RF_data_bn$Whitch_weed,
levels = c(0,1),
labels = c("Absent","Present"))
# Convert categorical variables to factors
RF_data_bn$IQR_soil_texture <- as.factor(RF_data_bn$IQR_soil_texture)
RF_data_bn$IQR_soil_color <- as.factor(RF_data_bn$IQR_soil_color)
RF_data_bn$IQF_agzone <- as.factor(RF_data_bn$IQF_agzone)
RF_data_bn$IQF_environment <- as.factor(RF_data_bn$IQF_environment)
RF_data_bn$IQF_position_slope <- as.factor(RF_data_bn$IQF_position_slope)
# Optional (often categorical as well)
RF_data_bn$IQR_plot_fert_quality <- as.factor(RF_data_bn$IQR_plot_fert_quality)Select variables and prepare dataset
library(tidyverse)
library(caret)
library(randomForest)
library(usdm)
# Select response + explanatory variables
data_model <- RF_data_bn %>%
dplyr::select(
avg_y_ratio_bn,
ICR_N_perc_23A,
#ICR_K_perc_23A,
K_factor,
Couch_Consit,
#Whitch_weed,
ICF_Alt_avg,
#ICR_Avail_P_avg,
#acc_rain_DAP120_mz,
#acc_rain_DAP30_bn,
acc_rain_DAP60_bn,
#acc_rain_DAP60_120_mz,
P_factor,
#ICR_arable_land_avg,
Comp_amount_corr_quali_bn,
Planting_date_rel_bn,
Slope,
#IQR_TF_tubura_client,
IQR_plot_fert_quality,
IQR_soil_texture,
ICR_Org_C_avg,
#IQR_soil_color,
IQF_position_slope,
#Weed_species_A_combined,
Weed_Overall_avg,
Mulch_planting_bn,
#IQF_agzone,
#IQF_environment
) %>%
na.omit()Train/ Test split
set.seed(42)
parts <- createDataPartition(data_model$avg_y_ratio_bn, p = 0.7, list = FALSE)
data_train <- data_model[parts, ]
data_test <- data_model[-parts, ]
x_train <- data_train[, -1]
y_train <- data_train$avg_y_ratio_bnVariable selection (RFE)
set.seed(627)
control_RFE <- rfeControl(
functions = rfFuncs,
method = "repeatedcv",
repeats = 1,
number = 5
)
start <- Sys.time()
result_RFE <- rfe(
x = x_train,
y = y_train,
sizes = 1:ncol(x_train),
rfeControl = control_RFE
)
end <- Sys.time()
print(end - start)## Time difference of 18.61518 secsSelected:
predictors(result_RFE)## [1] "Mulch_planting_bn" "acc_rain_DAP60_bn"
## [3] "ICF_Alt_avg" "Slope"
## [5] "Weed_Overall_avg" "Planting_date_rel_bn"
## [7] "IQR_soil_texture" "P_factor"
## [9] "ICR_Org_C_avg" "Comp_amount_corr_quali_bn"non selected
names(x_train)[!names(x_train) %in% result_RFE$optVariables]## [1] "ICR_N_perc_23A" "K_factor" "Couch_Consit"
## [4] "IQR_plot_fert_quality" "IQF_position_slope"Keep only selected variables
dataRFE <- subset(data_train, select = result_RFE$optVariables)I will force it to keep “Slope”, “IQR_soil_texture”, “acc_rain_p60_maize”, “ICR_Avail_P_avg”, “Comp_amount_corr_quali”, “IQF_position_slope”, “IQR_soil_texture”, “IQR_soil_color”
forced_vars <- c("Couch_Consit" )
final_vars <- unique(c(result_RFE$optVariables, forced_vars))
dataRFE <- data_train[, final_vars]Test colinearity
First we separate numeric from factoria as VIF only handles numeric
library(dplyr)
data_numeric <- dataRFE %>% dplyr::select(where(is.numeric))
data_factor <- dataRFE %>% dplyr::select(where(is.factor))
names(data_numeric)## [1] "Mulch_planting_bn" "acc_rain_DAP60_bn"
## [3] "ICF_Alt_avg" "Slope"
## [5] "Weed_Overall_avg" "Planting_date_rel_bn"
## [7] "ICR_Org_C_avg" "Comp_amount_corr_quali_bn"names(data_factor)## [1] "IQR_soil_texture" "Couch_Consit"Colineary test: VIF
Common thresholds:
VIF < 5 → safe VIF < 10 → acceptable VIF > 10 → problematic collinearity
So we keep all the numeric
library(usdm)
data_numeric_df <- as.data.frame(data_numeric)
vif_step <- vifstep(data_numeric_df, th = 10)
vif_step## No variable from the 8 input variables has collinearity problem.
##
## The linear correlation coefficients ranges between:
## min correlation ( Comp_amount_corr_quali_bn ~ Weed_Overall_avg ): -0.0007741101
## max correlation ( Planting_date_rel_bn ~ Mulch_planting_bn ): -0.4882837
##
## ---------- VIFs of the remained variables --------
## Variables VIF
## 1 Mulch_planting_bn 1.692869
## 2 acc_rain_DAP60_bn 1.215859
## 3 ICF_Alt_avg 1.951006
## 4 Slope 1.462851
## 5 Weed_Overall_avg 1.184360
## 6 Planting_date_rel_bn 1.459147
## 7 ICR_Org_C_avg 1.254915
## 8 Comp_amount_corr_quali_bn 1.109310Final RF database
data_train_final <- cbind(
avg_y_ratio = data_train$avg_y_ratio_bn,
dataRFE
)RF model calibration
# model calibration -------------------------------------------------------
fitControl <- trainControl(
method = "repeatedcv",
repeats = 1,
number = 5,
search = "grid"
)
# grid for ranger RF
rfGrid <- expand.grid(
mtry = c(1:length(dataRFE)),
min.node.size = c(5, 10, 50, 100),
splitrule = "variance"
)
set.seed(627)
start <- Sys.time()
RF_model_train <- train(
formula(paste0("avg_y_ratio_bn ~ ",
paste0(names(dataRFE), collapse = " + "))),
data = data_train,
method = "ranger",
trControl = fitControl,
tuneGrid = rfGrid,
metric = "RMSE",
verbose = FALSE
)
end <- Sys.time()
print(end - start)## Time difference of 10.2635 secsFinal model using best hyperparameters
fitControl <- trainControl(method = "none", search = "grid")
RFFGrid <- expand.grid(
mtry = RF_model_train$bestTune$mtry,
splitrule = RF_model_train$bestTune$splitrule,
min.node.size = RF_model_train$bestTune$min.node.size
)
RFFGrid## mtry splitrule min.node.size
## 1 2 variance 5Train final RF
fit_RF <- train(
formula(paste0("avg_y_ratio_bn ~ ",
paste0(names(dataRFE), collapse = " + "))),
data = data_train,
method = "ranger",
trControl = fitControl,
metric = "RMSE",
verbose = FALSE,
tuneGrid = RFFGrid
)Train set
library(Metrics)
results_RF_train <- data.frame(
avg_y_ratio_bn = data_train$avg_y_ratio_bn,
y_pred = predict(fit_RF, newdata = data_train)
)
results_RF_train$residual <-
results_RF_train$avg_y_ratio_bn - results_RF_train$y_pred
RMSE_RF_train <- Metrics::rmse(results_RF_train$avg_y_ratio_bn, results_RF_train$y_pred)
RMSE_RF_train## [1] 0.1047524R2_RF_train <- cor(results_RF_train$y_pred, results_RF_train$avg_y_ratio_bn)^2
R2_RF_train## [1] 0.9036525For the test set:
results_RF_test <- data.frame(
avg_y_ratio_bn = data_test$avg_y_ratio_bn,
y_pred = predict(fit_RF, newdata = data_test)
)
results_RF_test$residual <- results_RF_test$avg_y_ratio_bn - results_RF_test$y_pred
RMSE_RF_test <- Metrics::rmse(results_RF_test$avg_y_ratio_bn, results_RF_test$y_pred)
RMSE_RF_test## [1] 0.1849755R2_RF_test <- cor(results_RF_test$y_pred, results_RF_test$avg_y_ratio_bn)^2
R2_RF_test## [1] 0.2706352Observed vs Predicted plots
Trained data
plot(results_RF_train$y_pred, results_RF_train$avg_y_ratio_bn,
xlab = "Predicted Y ratio",
ylab = "Observed Y ratio",
pch = 20,
main = "RF Train")
abline(lm(avg_y_ratio_bn ~ y_pred, data = results_RF_train), col = "red")
abline(0,1,col="red", lty=2)
legend("topleft",
legend = c(paste("R2 =", round(R2_RF_train,3)),
paste("RMSE =", round(RMSE_RF_train,3))),
bty="n")
tested
# Counts
true_good <- sum(results_RF_test$avg_y_ratio_bn >= 0.9 & results_RF_test$y_pred >= 0.9)
false_acc <- sum(results_RF_test$avg_y_ratio_bn < 0.9 & results_RF_test$y_pred >= 0.9)
plot(results_RF_test$y_pred, results_RF_test$avg_y_ratio_bn,
xlab = "Predicted CA:CT yield ratio",
ylab = "Observed CA:CT yield ratio",
pch = 20,
main = "RF Test")
abline(lm(avg_y_ratio_bn ~ y_pred, data = results_RF_test), col = "red")
abline(0,1,col="red", lty=2)
# Threshold lines
abline(v = 0.9, col = "blue", lty = 2)
abline(h = 0.9, col = "blue", lty = 2)
# Add quadrant labels
text(x = 0.8, y = 1.2,
labels = paste("True performer=", true_good),
col = "darkgreen", pos = 4)
text(x = 0.8, y = 0.5,
labels = paste("False performer=", false_acc),
col = "orange", pos = 4)
legend("topleft",
legend = c(paste("R² =", round(R2_RF_test,3)),
paste("RMSE =", round(RMSE_RF_test,3))),
bty="n")
Stave model objects
save(
data_model,
parts, data_train, data_test,
result_RFE,
RF_model_train, fit_RF,
results_RF_train, RMSE_RF_train, R2_RF_train,
results_RF_test, RMSE_RF_test, R2_RF_test,
file = "RF_avg_Y_ratio_model.Rdata"
)names(dataRFE) %in% names(data_train)
names(dataRFE) %in% names(data_train)## [1] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUEVisualizetion of variables importance
# =========================================
# Variable Importance (Beans - FINAL)
# =========================================
library(dplyr)
library(stringr)
library(ggplot2)
library(caret)
# -----------------------------
# 1. Train model WITH importance
# -----------------------------
fit_RF <- train(
formula(paste0("avg_y_ratio_bn ~ ",
paste0(names(dataRFE), collapse = " + "))),
data = data_train,
method = "ranger",
trControl = trainControl(method = "none"),
tuneGrid = RFFGrid,
importance = "permutation" # 🔥 CRITICAL
)
# -----------------------------
# 2. Extract importance
# -----------------------------
importance_RF <- varImp(fit_RF)
imp <- importance_RF$importance
imp$variable <- rownames(imp)
# -----------------------------
# 3. Aggregate categories → variables
# -----------------------------
imp_grouped <- imp %>%
mutate(
base_var = case_when(
str_detect(variable, "^K_factor") ~ "K_factor",
str_detect(variable, "^P_factor") ~ "P_factor",
str_detect(variable, "^Couch_Consit") ~ "Couch_Consit",
str_detect(variable, "^Whitch_weed") ~ "Whitch_weed",
str_detect(variable, "^IQR_soil_texture") ~ "IQR_soil_texture",
str_detect(variable, "^IQR_soil_color") ~ "IQR_soil_color",
str_detect(variable, "^IQF_position_slope") ~ "IQF_position_slope",
str_detect(variable, "^IQR_plot_fert_quality") ~ "IQR_plot_fert_quality",
str_detect(variable, "^Weed_species_A_combined") ~ "Weed_species_A_combined",
TRUE ~ variable
)
) %>%
group_by(base_var) %>%
summarise(
Overall = sum(Overall),
.groups = "drop"
) %>%
arrange(desc(Overall))
# -----------------------------
# 4. Plot (same style as yours)
# -----------------------------
ggplot(imp_grouped, aes(x = reorder(base_var, Overall), y = Overall)) +
geom_col() +
coord_flip() +
theme_minimal() +
labs(
title = "Variable Importance (Beans - Aggregated)",
x = "Variable",
y = "Overall"
)