— Databases
# Load required package for Excel files
library(readxl)
# Define the raw data folder path
data_raw <- file.path(getwd(), "data_raw")
# --- File 1: X26B_CAs_farmers.xlsx ---
selection_CA1 <- read_excel(
file.path(data_raw, "26B CAs farmers.xlsx"),
sheet = "26B CA2_Mz | Bean rotation"
)
## New names:
## • `Region` -> `Region...2`
## • `` -> `...14`
## • `` -> `...15`
## • `` -> `...16`
## • `Category` -> `Category...17`
## • `Nbr of farmers` -> `Nbr of farmers...18`
## • `` -> `...19`
## • `Region` -> `Region...20`
## • `Category` -> `Category...21`
## • `Nbr of farmers` -> `Nbr of farmers...22`
selection_CA2 <- read_excel(
file.path(data_raw, "26B CAs farmers.xlsx"),
sheet = "26B CA3_MzBean IC"
)
## New names:
## • `` -> `...14`
## • `` -> `...17`
# --- File 2: CONV_CA_Ratio.csv ---
yield_CA1 <- read.csv(
file.path(data_raw, "CONV_CA_Ratio.csv"),
stringsAsFactors = FALSE
)
# --- File 3: CA3 data.xlsx (3 sheets) ---
yield_CA2_25A <- read_excel(
file.path(data_raw, "CA3 data.xlsx"),
sheet = "25A CA3"
)
## New names:
## • `seed grams planted` -> `seed grams planted...29`
## • `seed grams planted` -> `seed grams planted...31`
## • `binnage2done` -> `binnage2done...149`
## • `binnage2done` -> `binnage2done...150`
## • `rustpresent.x` -> `rustpresent.x...235`
## • `rustpresent.x` -> `rustpresent.x...245`
## • `numberplantssampletaken` -> `numberplantssampletaken...247`
## • `numberplantssampletaken` -> `numberplantssampletaken...248`
## • `weedsreduceharvest` -> `weedsreduceharvest...293`
## • `weedfirstspecies` -> `weedfirstspecies...294`
## • `othermainspecies` -> `othermainspecies...295`
## • `firstspeciespicture` -> `firstspeciespicture...296`
## • `firstspeciespictureURL` -> `firstspeciespictureURL...297`
## • `weedsecondspecies` -> `weedsecondspecies...298`
## • `secondspeciesother` -> `secondspeciesother...299`
## • `secondspeciespicture` -> `secondspeciespicture...300`
## • `secondspeciespictureURL` -> `secondspeciespictureURL...301`
## • `greencoverpercentageharvest` -> `greencoverpercentageharvest...302`
## • `greencoverpicture` -> `greencoverpicture...303`
## • `greencoverpictureURL` -> `greencoverpictureURL...304`
## • `weedsreduceharvest` -> `weedsreduceharvest...306`
## • `weedfirstspecies` -> `weedfirstspecies...307`
## • `othermainspecies` -> `othermainspecies...308`
## • `firstspeciespicture` -> `firstspeciespicture...309`
## • `firstspeciespictureURL` -> `firstspeciespictureURL...310`
## • `weedsecondspecies` -> `weedsecondspecies...311`
## • `secondspeciesother` -> `secondspeciesother...312`
## • `secondspeciespicture` -> `secondspeciespicture...313`
## • `secondspeciespictureURL` -> `secondspeciespictureURL...314`
## • `greencoverpercentageharvest` -> `greencoverpercentageharvest...315`
## • `greencoverpicture` -> `greencoverpicture...316`
## • `greencoverpictureURL` -> `greencoverpictureURL...317`
## • `questionswantmention` -> `questionswantmention...323`
## • `howmuchreduceyield` -> `howmuchreduceyield...325`
## • `maize_ green_ biomass_ total_ weight_ at_ harvest` -> `maize_ green_
## biomass_ total_ weight_ at_ harvest...385`
## • `maize_ green_ biomass_ total_ weight_ at_ harvest` -> `maize_ green_
## biomass_ total_ weight_ at_ harvest...388`
## • `questionswantmention` -> `questionswantmention...407`
## • `howmuchreduceyield` -> `howmuchreduceyield...408`
yield_CA2_25B <- read_excel(
file.path(data_raw, "CA3 data.xlsx"),
sheet = " 25B CA3"
)
## New names:
## • `environment` -> `environment...3`
## • `Spacing_btwn_rows_cm` -> `Spacing_btwn_rows_cm...31`
## • `Spacing_btwn_rows_cm` -> `Spacing_btwn_rows_cm...32`
## • `spacing_between_holes` -> `spacing_between_holes...33`
## • `spacing_between_holes` -> `spacing_between_holes...34`
## • `land_prep_sex` -> `land_prep_sex...41`
## • `land_prep_sex` -> `land_prep_sex...54`
## • `urea_application_period` -> `urea_application_period...76`
## • `topdress_grams_used` -> `topdress_grams_used...77`
## • `topdress_application_date` -> `topdress_application_date...78`
## • `topdress_green_cover_percentage` -> `topdress_green_cover_percentage...79`
## • `labor_start_topdress` -> `labor_start_topdress...80`
## • `labor_end_topdress` -> `labor_end_topdress...81`
## • `topdress_minutes` -> `topdress_minutes...82`
## • `labor_start_binnage` -> `labor_start_binnage...90`
## • `labor_end_binnage` -> `labor_end_binnage...91`
## • `binnage_minutes` -> `binnage_minutes...92`
## • `binnage_gender` -> `binnage_gender...93`
## • `binnage_easy` -> `binnage_easy...95`
## • `binnage_easy_compared_CT` -> `binnage_easy_compared_CT...96`
## • `binnage_recommandation` -> `binnage_recommandation...97`
## • `binnage_recommandation_explaination` ->
## `binnage_recommandation_explaination...98`
## • `binnage_recommandation_marks` -> `binnage_recommandation_marks...99`
## • `binnage_recommandation_marks_explaination` ->
## `binnage_recommandation_marks_explaination...100`
## • `binnage_complication_observed` -> `binnage_complication_observed...101`
## • `Alternative_weedrec_ethod_by_farmer` ->
## `Alternative_weedrec_ethod_by_farmer...102`
## • `binnage_done_differently` -> `binnage_done_differently...103`
## • `binnage2_done...401` -> `binnage2_done`
## • `binnage2_date_...405` -> `binnage2_date_...105`
## • `binnage2_io_present_...409` -> `binnage2_io_present_...106`
## • `binnage2_method_...413` -> `binnage2_method_...107`
## • `labor_start_binnage` -> `labor_start_binnage...108`
## • `labor_end_binnage` -> `labor_end_binnage...109`
## • `binnage_minutes` -> `binnage_minutes...110`
## • `binnage_gender` -> `binnage_gender...111`
## • `binnage_easy` -> `binnage_easy...113`
## • `binnage_easy_compared_CT` -> `binnage_easy_compared_CT...114`
## • `binnage_recommandation` -> `binnage_recommandation...115`
## • `binnage_recommandation_explaination` ->
## `binnage_recommandation_explaination...116`
## • `binnage_recommandation_marks` -> `binnage_recommandation_marks...117`
## • `binnage_recommandation_marks_explaination` ->
## `binnage_recommandation_marks_explaination...118`
## • `binnage_complication_observed` -> `binnage_complication_observed...119`
## • `binnage_done_differently` -> `binnage_done_differently...120`
## • `environment` -> `environment...121`
## • `binnage2_done_...569` -> `binnage2_done_`
## • `binnage2_date_...573` -> `binnage2_date_...123`
## • `binnage2_io_present_...577` -> `binnage2_io_present_...124`
## • `binnage2_method_...581` -> `binnage2_method_...125`
## • `Alternative_weedrec_ethod_by_farmer` ->
## `Alternative_weedrec_ethod_by_farmer...126`
## • `binnage3_done_...590` -> `binnage3_done_...127`
## • `binnage3_date_...594` -> `binnage3_date_`
## • `binnage3_method_...602` -> `binnage3_method_...130`
## • `labor_start_binnage` -> `labor_start_binnage...131`
## • `labor_end_binnage` -> `labor_end_binnage...132`
## • `binnage_minutes` -> `binnage_minutes...133`
## • `binnage_gender` -> `binnage_gender...134`
## • `binnage_easy` -> `binnage_easy...136`
## • `binnage_easy_compared_CT` -> `binnage_easy_compared_CT...137`
## • `binnage_recommandation` -> `binnage_recommandation...138`
## • `binnage_recommandation_explaination` ->
## `binnage_recommandation_explaination...139`
## • `binnage_recommandation_marks` -> `binnage_recommandation_marks...140`
## • `binnage_recommandation_marks_explaination` ->
## `binnage_recommandation_marks_explaination...141`
## • `binnage_complication_observed` -> `binnage_complication_observed...142`
## • `binnage_done_differently` -> `binnage_done_differently...143`
## • `binnage3_done_...758` -> `binnage3_done_...145`
## • `binnage3_io_present_...766` -> `binnage3_io_present_`
## • `binnage3_method_...770` -> `binnage3_method_...148`
## • `Alternative_weedrec_ethod_by_farmer` ->
## `Alternative_weedrec_ethod_by_farmer...149`
## • `urea_application_period` -> `urea_application_period...150`
## • `topdress_grams_used` -> `topdress_grams_used...151`
## • `topdress_application_date` -> `topdress_application_date...152`
## • `topdress_green_cover_percentage` -> `topdress_green_cover_percentage...153`
## • `labor_start_topdress` -> `labor_start_topdress...154`
## • `labor_end_topdress` -> `labor_end_topdress...155`
## • `topdress_minutes` -> `topdress_minutes...156`
## • `date_recording_disease` -> `date_recording_disease...166`
## • `beans_anthracnose_incidence` -> `beans_anthracnose_incidence...167`
## • `beans_anthracnose_severity` -> `beans_anthracnose_severity...168`
## • `beans_anthracnose_plantcount` -> `beans_anthracnose_plantcount...169`
## • `angular_leaf_spot_present` -> `angular_leaf_spot_present...170`
## • `desease_name_angular_leaf_spot` -> `desease_name_angular_leaf_spot...171`
## • `angular_leaf_spot_incidence` -> `angular_leaf_spot_incidence...173`
## • `angular_leaf_spot_severity` -> `angular_leaf_spot_severity...174`
## • `angular_leaf_spot_plantcount` -> `angular_leaf_spot_plantcount...175`
## • `beans_anthracnose_incidence` -> `beans_anthracnose_incidence...199`
## • `beans_anthracnose_severity` -> `beans_anthracnose_severity...200`
## • `beans_anthracnose_plantcount` -> `beans_anthracnose_plantcount...201`
## • `angular_leaf_spot_present` -> `angular_leaf_spot_present...202`
## • `desease_name_angular_leaf_spot` -> `desease_name_angular_leaf_spot...203`
## • `angular_leaf_spot_incidence` -> `angular_leaf_spot_incidence...204`
## • `angular_leaf_spot_severity` -> `angular_leaf_spot_severity...205`
## • `angular_leaf_spot_plantcount` -> `angular_leaf_spot_plantcount...206`
## • `date_recording_disease` -> `date_recording_disease...207`
## • `rust_recording-date` -> `rust_recording-date...240`
## • `gls_severity` -> `gls_severity...241`
## • `gls_recording-date` -> `gls_recording-date...242`
## • `rust_present` -> `rust_present...243`
## • `rust_Incidence` -> `rust_Incidence...244`
## • `rust_severity` -> `rust_severity...245`
## • `tlb_recording-date` -> `tlb_recording-date...249`
## • `gls_severity` -> `gls_severity...252`
## • `rust_present` -> `rust_present...253`
## • `rust_Incidence` -> `rust_Incidence...254`
## • `rust_severity` -> `rust_severity...255`
## • `rust_recording-date` -> `rust_recording-date...256`
## • `gls_recording-date` -> `gls_recording-date...257`
## • `tlb_recording-date` -> `tlb_recording-date...258`
## • `living_much_day` -> `living_much_day...259`
## • `weeds_reduce_harvest` -> `weeds_reduce_harvest...260`
## • `weed_first_species` -> `weed_first_species...261`
## • `other_main_species` -> `other_main_species...262`
## • `first_species_picture` -> `first_species_picture...263`
## • `weed_second_species` -> `weed_second_species...265`
## • `second_species_other` -> `second_species_other...266`
## • `second_species_picture` -> `second_species_picture...267`
## • `green_cover_percentage` -> `green_cover_percentage...269`
## • `Green_cover_picture` -> `Green_cover_picture...270`
## • `living_much_day` -> `living_much_day...272`
## • `weeds_reduce_harvest` -> `weeds_reduce_harvest...273`
## • `weed_first_species` -> `weed_first_species...274`
## • `other_main_species` -> `other_main_species...275`
## • `first_species_picture` -> `first_species_picture...276`
## • `weed_second_species` -> `weed_second_species...278`
## • `second_species_other` -> `second_species_other...279`
## • `second_species_picture` -> `second_species_picture...280`
## • `green_cover_percentage` -> `green_cover_percentage...282`
## • `Green_cover_picture` -> `Green_cover_picture...283`
## • `questions_want_mention` -> `questions_want_mention...287`
## • `questions_affect_control` -> `questions_affect_control...288`
## • `Officer_do_you_thin_affect_control` ->
## `Officer_do_you_thin_affect_control...289`
## • `How_much_it_will_reduce_yield` -> `How_much_it_will_reduce_yield...290`
## • `Harvest_box_length_m` -> `Harvest_box_length_m...293`
## • `Harvest_box_width_m` -> `Harvest_box_width_m...294`
## • `Harvest_box_size_m2` -> `Harvest_box_size_m2...295`
## • `residues_mulch_reduced_season` -> `residues_mulch_reduced_season...304`
## • `residues_mulch_reduced_reason_harvest` ->
## `residues_mulch_reduced_reason_harvest...305`
## • `others_reasons_ask_farm` -> `others_reasons_ask_farm...306`
## • `questions_want_mention` -> `questions_want_mention...307`
## • `questions_affect_control` -> `questions_affect_control...308`
## • `how_much_reduce_yield` -> `how_much_reduce_yield...309`
## • `revenue_possible_harvest_box` -> `revenue_possible_harvest_box...311`
## • `beans_green_biomass_total_weig` -> `beans_green_biomass_total_weig...314`
## • `beans_green_biomass_sample_weight` ->
## `beans_green_biomass_sample_weight...318`
## • `beans_green_biomass_sample_code` -> `beans_green_biomass_sample_code...319`
## • `harvest_methods_recommadation` -> `harvest_methods_recommadation...320`
## • `harvest_method_recomendation_reasons` ->
## `harvest_method_recomendation_reasons...321`
## • `Moisture_content` -> `Moisture_content...324`
## • `residues_mulch_reduced_season` -> `residues_mulch_reduced_season...325`
## • `residues_mulch_reduced_reason_harvest` ->
## `residues_mulch_reduced_reason_harvest...326`
## • `others_reasons_ask_farm` -> `others_reasons_ask_farm...327`
## • `questions_want_mention` -> `questions_want_mention...328`
## • `questions_affect_control` -> `questions_affect_control...329`
## • `how_much_reduce_yield` -> `how_much_reduce_yield...330`
## • `revenue_possible_harvest_box` -> `revenue_possible_harvest_box...331`
## • `beans_green_biomass_total_weig` -> `beans_green_biomass_total_weig...332`
## • `beans_green_biomass_total_weight` ->
## `beans_green_biomass_total_weight...334`
## • `beans_green_biomass_total_weight` ->
## `beans_green_biomass_total_weight...336`
## • `beans_green_biomass_sample_weight` ->
## `beans_green_biomass_sample_weight...338`
## • `beans_green_biomass_sample_code` -> `beans_green_biomass_sample_code...339`
## • `harvest_methods_recommadation` -> `harvest_methods_recommadation...340`
## • `harvest_method_recomendation_reasons` ->
## `harvest_method_recomendation_reasons...341`
## • `Officer_do_you_thin_affect_control` ->
## `Officer_do_you_thin_affect_control...372`
## • `How_much_it_will_reduce_yield` -> `How_much_it_will_reduce_yield...373`
## • `Harvest_box_length_m` -> `Harvest_box_length_m...375`
## • `Harvest_box_width_m` -> `Harvest_box_width_m...376`
## • `Harvest_box_size_m2` -> `Harvest_box_size_m2...377`
## • `Moisture_content` -> `Moisture_content...400`
yield_CA2_26A <- read_excel(
file.path(data_raw, "CA3 data.xlsx"),
sheet = "26A CA3"
)
## New names:
## • `q085b_prob1_happened...340` -> `q085b_prob1_happened...336`
## • `q084d_prob2_yield_reduction...351` -> `q084d_prob2_yield_reduction...347`
## • `q085b_prob1_happened...549` -> `q085b_prob1_happened...525`
## • `q084d_prob2_yield_reduction...551` -> `q084d_prob2_yield_reduction...527`
— Balance check in old CA trial
We firs merge databases
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, union
# --- 1. Confirm column names exist ---
cols_wanted_sel <- c("District", "Agzone", "Officer", "Cell", "Block", "Quali")
cols_wanted_yield <- c("IQR_Season", "IQF_environment", "IQR_sex_HH_head", "Block",
"ICF_Lat", "ICF_Long", "ICT_Alt", "Y_ratio", "env_index",
"yield_A", "yield_C", "IQF_position_slope", "Erosion_quanti", "Slope")
cat("Missing in selection_CA1:\n")
## Missing in selection_CA1:
print(cols_wanted_sel[!cols_wanted_sel %in% names(selection_CA1)])
## [1] "Cell"
cat("\nMissing in yield_CA1:\n")
##
## Missing in yield_CA1:
print(cols_wanted_yield[!cols_wanted_yield %in% names(yield_CA1)])
## [1] "ICT_Alt"
# --- 2. Compare Block values between the two datasets ---
blocks_sel <- unique(selection_CA1$Block)
blocks_yield <- unique(yield_CA1$Block)
cat("\nNumber of unique Blocks in selection_CA1:", length(blocks_sel), "\n")
##
## Number of unique Blocks in selection_CA1: 266
cat("Number of unique Blocks in yield_CA1:", length(blocks_yield), "\n")
## Number of unique Blocks in yield_CA1: 565
# Blocks listed in selection_CA1 but NOT found anywhere in yield_CA1 --> PROBLEM
orphan_blocks <- setdiff(blocks_sel, blocks_yield)
cat("\nNumber of Blocks in selection_CA1 with NO match in yield_CA1:", length(orphan_blocks), "\n")
##
## Number of Blocks in selection_CA1 with NO match in yield_CA1: 3
if (length(orphan_blocks) > 0) {
cat("These Blocks are missing from yield_CA1:\n")
print(orphan_blocks)
}
## These Blocks are missing from yield_CA1:
## [1] "48_10" "115_4" "96_16"
# Optional: view the actual selection_CA1 rows that are orphaned, for context
orphan_rows <- selection_CA1 %>% filter(Block %in% orphan_blocks)
print(orphan_rows)
## # A tibble: 3 × 22
## Season Region...2 Agzone District Officer cell `Trial number` farmer_id Block
## <chr> <chr> <chr> <chr> <chr> <chr> <dbl> <dbl> <chr>
## 1 26B NW Volca… Burera Uwimana Nkum… 48 10 48_10
## 2 26B SW Cyang… Rusizi Vedaste Gahi… 115 4 115_4
## 3 26B W Lake … Nyamash… Joseli… Jara… 96 16 96_16
## # ℹ 13 more variables: Farmer_name <chr>, village <chr>, phone <dbl>,
## # Quali <chr>, ...14 <lgl>, ...15 <lgl>, ...16 <lgl>, Category...17 <chr>,
## # `Nbr of farmers...18` <dbl>, ...19 <dbl>, Region...20 <chr>,
## # Category...21 <chr>, `Nbr of farmers...22` <dbl>
library(dplyr)
# --- Build check_CA1 ---
# Subset selection_CA1 (note: lowercase "cell")
sel_subset <- selection_CA1 %>%
select(District, Agzone, Officer, cell, Block, Quali)
# Subset yield_CA1 (note: ICF_Alt, not ICT_Alt)
yield_subset <- yield_CA1 %>%
select(Block, IQR_Season, IQF_environment, IQR_sex_HH_head,
ICF_Lat, ICF_Long, ICF_Alt, Y_ratio, env_index, yield_A, yield_C)
# Left join: keep ALL Blocks from selection_CA1, even if no match in yield_CA1
# This naturally keeps all seasons per Block where they exist
check_CA1 <- sel_subset %>%
left_join(yield_subset, by = "Block") %>%
mutate(
data_flag = ifelse(is.na(IQR_Season), "no_yield_data_found", "ok")
)
## Warning in left_join(., yield_subset, by = "Block"): Detected an unexpected many-to-many relationship between `x` and `y`.
## ℹ Row 1 of `x` matches multiple rows in `y`.
## ℹ Row 1855 of `y` matches multiple rows in `x`.
## ℹ If a many-to-many relationship is expected, set `relationship =
## "many-to-many"` to silence this warning.
# --- Sanity checks ---
cat("Total rows in check_CA1:", nrow(check_CA1), "\n")
## Total rows in check_CA1: 1527
cat("Unique Blocks in check_CA1:", length(unique(check_CA1$Block)), "\n")
## Unique Blocks in check_CA1: 266
cat("Rows flagged as missing yield data:", sum(check_CA1$data_flag == "no_yield_data_found"), "\n\n")
## Rows flagged as missing yield data: 3
# View the flagged/orphan rows specifically
check_CA1 %>% filter(data_flag == "no_yield_data_found")
## # A tibble: 3 × 17
## District Agzone Officer cell Block Quali IQR_Season IQF_environment
## <chr> <chr> <chr> <chr> <chr> <chr> <chr> <chr>
## 1 Burera Volcanic cones Uwimana Nkum… 48_10 Very… <NA> <NA>
## 2 Rusizi Cyangugu Vedaste Gahi… 115_4 Very… <NA> <NA>
## 3 Nyamasheke Lake Kivu Joseli… Jara… 96_16 Very… <NA> <NA>
## # ℹ 9 more variables: IQR_sex_HH_head <chr>, ICF_Lat <dbl>, ICF_Long <dbl>,
## # ICF_Alt <dbl>, Y_ratio <dbl>, env_index <dbl>, yield_A <dbl>,
## # yield_C <dbl>, data_flag <chr>
# How many seasons per Block did we retain (for the non-orphan Blocks)?
check_CA1 %>%
filter(data_flag == "ok") %>%
count(Block, name = "n_seasons") %>%
arrange(desc(n_seasons)) %>%
head(10)
## # A tibble: 10 × 2
## Block n_seasons
## <chr> <int>
## 1 16_15 12
## 2 22_2 12
## 3 17_2 8
## 4 102_10 6
## 5 102_11 6
## 6 102_12 6
## 7 102_13 6
## 8 102_14 6
## 9 102_2 6
## 10 102_5 6
Then explore the balce in terms of yield
library(dplyr)
library(tidyr)
library(ggplot2)
# --- Filter out orphan/no-data rows ---
check_CA1_clean <- check_CA1 %>%
filter(data_flag != "no_yield_data_found")
# --- Order Quali levels logically ---
check_CA1_clean <- check_CA1_clean %>%
mutate(Quali = factor(Quali, levels = c("Very Low", "Low", "Medium", "High", "Very High")))
# --- Summary table: mean Y_ratio, yield_A, yield_C by Quali x Season ---
summary_table <- check_CA1_clean %>%
group_by(Quali, IQR_Season) %>%
summarise(
n = n(),
mean_Y_ratio = mean(Y_ratio, na.rm = TRUE),
mean_yield_A = mean(yield_A, na.rm = TRUE),
mean_yield_C = mean(yield_C, na.rm = TRUE),
.groups = "drop"
) %>%
arrange(IQR_Season, Quali)
print(summary_table)
## # A tibble: 30 × 6
## Quali IQR_Season n mean_Y_ratio mean_yield_A mean_yield_C
## <fct> <chr> <int> <dbl> <dbl> <dbl>
## 1 Very Low 23A 65 0.932 4.44 3.98
## 2 Low 23A 11 0.840 5.59 4.57
## 3 Medium 23A 39 0.975 5.22 4.82
## 4 High 23A 5 0.913 5.39 4.93
## 5 Very High 23A 133 0.919 5.27 4.70
## 6 Very Low 23B 64 0.780 1.39 1.01
## 7 Low 23B 11 0.971 1.28 0.989
## 8 Medium 23B 41 0.867 1.38 1.18
## 9 High 23B 5 0.729 1.63 1.26
## 10 Very High 23B 139 0.763 1.75 1.28
## # ℹ 20 more rows
# Optional: save table to outputs folder
write.csv(summary_table, file.path(output_dir, "summary_Quali_Season.csv"), row.names = FALSE)
library(dplyr)
library(ggplot2)
library(multcompView)
# --- Function: run ANOVA + Tukey HSD per season, return letters ---
get_letters <- function(df) {
model <- aov(Y_ratio ~ Quali, data = df)
tukey <- TukeyHSD(model)
letters_df <- multcompLetters4(model, tukey)$Quali$Letters
data.frame(Quali = names(letters_df), letter = letters_df, row.names = NULL)
}
# --- Letters per season ---
letters_by_season <- check_CA1_clean %>%
filter(!is.na(Y_ratio)) %>%
group_by(IQR_Season) %>%
group_modify(~ get_letters(.x)) %>%
ungroup()
# --- N, mean, and y-position per Quali x Season ---
label_stats <- check_CA1_clean %>%
filter(!is.na(Y_ratio)) %>%
group_by(IQR_Season, Quali) %>%
summarise(
n = n(),
mean_v = mean(Y_ratio, na.rm = TRUE),
y_max = max(Y_ratio, na.rm = TRUE),
.groups = "drop"
)
# --- Combine letters + stats into one label, stacked with line breaks ---
label_data <- label_stats %>%
left_join(letters_by_season, by = c("IQR_Season", "Quali")) %>%
mutate(
label = paste0(letter, "\nn=", n, "\nmean=", round(mean_v, 2)),
y_pos = y_max * 1.05
)
# --- Final plot ---
p_yratio <- ggplot(check_CA1_clean %>% filter(!is.na(Y_ratio)),
aes(x = Quali, y = Y_ratio, fill = Quali)) +
geom_boxplot(na.rm = TRUE) +
facet_wrap(~ IQR_Season) +
geom_text(data = label_data,
aes(x = Quali, y = y_pos, label = label),
inherit.aes = FALSE, vjust = 0, size = 3, fontface = "bold", lineheight = 0.85) +
labs(title = "Y_ratio by Quali, faceted by Season (letters = Tukey HSD groups)",
x = "Quali", y = "Y_ratio") +
theme_minimal() +
theme(axis.text.x = element_text(angle = 45, hjust = 1), legend.position = "none") +
scale_y_continuous(expand = expansion(mult = c(0.05, 0.2))) # extra headroom for labels
print(p_yratio)
ggsave(file.path(output_dir, "boxplot_Y_ratio_by_Quali_Season_letters.png"), p_yratio, width = 9, height = 6)
across seasons
library(dplyr)
library(ggplot2)
library(multcompView)
# --- Step 1: average Y_ratio per Block across seasons ---
block_avg <- check_CA1_clean %>%
group_by(Block, District, Agzone, Officer, cell, Quali) %>%
summarise(
n_seasons = n(),
mean_Y_ratio = mean(Y_ratio, na.rm = TRUE),
.groups = "drop"
)
cat("Number of Blocks after averaging:", nrow(block_avg), "\n")
## Number of Blocks after averaging: 264
head(block_avg)
## # A tibble: 6 × 8
## Block District Agzone Officer cell Quali n_seasons mean_Y_ratio
## <chr> <chr> <chr> <chr> <chr> <fct> <int> <dbl>
## 1 102_10 Bugarama Bugarama Wellars Hangabashi Very Low 6 0.647
## 2 102_11 Bugarama Bugarama Wellars Hangabashi Very Low 6 0.996
## 3 102_12 Bugarama Bugarama Wellars Hangabashi Very High 6 0.658
## 4 102_13 Bugarama Bugarama Wellars Hangabashi Very Low 6 0.890
## 5 102_14 Bugarama Bugarama Wellars Hangabashi Very High 6 1.12
## 6 102_2 Bugarama Bugarama Wellars Hangabashi Very High 6 0.764
# --- Step 2: Tukey HSD across Quali groups (no season facet now) ---
model <- aov(mean_Y_ratio ~ Quali, data = block_avg)
tukey <- TukeyHSD(model)
letters_df <- multcompLetters4(model, tukey)$Quali$Letters
letters_tbl <- data.frame(Quali = names(letters_df), letter = letters_df, row.names = NULL)
# --- Step 3: N, mean, y-position per Quali group ---
label_stats <- block_avg %>%
group_by(Quali) %>%
summarise(
n = n(),
mean_v = mean(mean_Y_ratio, na.rm = TRUE),
y_max = max(mean_Y_ratio, na.rm = TRUE),
.groups = "drop"
)
label_data <- label_stats %>%
left_join(letters_tbl, by = "Quali") %>%
mutate(
label = paste0(letter, "\nn=", n, "\nmean=", round(mean_v, 2)),
y_pos = y_max * 1.05
)
# --- Step 4: Plot ---
library(dplyr)
library(ggplot2)
library(multcompView)
p_block_avg <- ggplot(block_avg, aes(x = Quali, y = mean_Y_ratio, fill = Quali)) +
geom_boxplot(na.rm = TRUE, outlier.shape = NA, alpha = 0.6) + # outlier.shape = NA avoids double-plotting outliers
geom_jitter(width = 0.15, height = 0, alpha = 0.6, size = 1.8, color = "black") +
geom_text(data = label_data,
aes(x = Quali, y = y_pos, label = label),
inherit.aes = FALSE, vjust = 0, size = 3, fontface = "bold", lineheight = 0.85) +
labs(title = "Average Y_ratio per Block (across seasons) by Quali",
x = "Quality of the trial", y = "Mean Y_ratio (per Block)") +
theme_minimal() +
theme(axis.text.x = element_text(angle = 45, hjust = 1), legend.position = "none") +
scale_y_continuous(expand = expansion(mult = c(0.05, 0.2)))
print(p_block_avg)
ggsave(file.path(output_dir, "boxplot_Y_ratio_by_Block_avg_Quali_dots.png"), p_block_avg, width = 8, height = 6)
Effect of selection of slope
library(dplyr)
library(ggplot2)
library(multcompView)
# --- Step 0: rebuild check_CA1 to include Slope ---
yield_subset <- yield_CA1 %>%
select(Block, IQR_Season, IQF_environment, IQR_sex_HH_head,
ICF_Lat, ICF_Long, ICF_Alt, Y_ratio, env_index, yield_A, yield_C, Slope)
check_CA1 <- sel_subset %>%
left_join(yield_subset, by = "Block") %>%
mutate(
data_flag = ifelse(is.na(IQR_Season), "no_yield_data_found", "ok")
)
## Warning in left_join(., yield_subset, by = "Block"): Detected an unexpected many-to-many relationship between `x` and `y`.
## ℹ Row 1 of `x` matches multiple rows in `y`.
## ℹ Row 1855 of `y` matches multiple rows in `x`.
## ℹ If a many-to-many relationship is expected, set `relationship =
## "many-to-many"` to silence this warning.
check_CA1_clean <- check_CA1 %>%
filter(data_flag != "no_yield_data_found") %>%
mutate(Quali = factor(Quali, levels = c("Very Low", "Low", "Medium", "High", "Very High")))
# --- Step 1: average Slope per Block across seasons ---
block_avg_slope <- check_CA1_clean %>%
group_by(Block, District, Agzone, Officer, cell, Quali) %>%
summarise(
n_seasons = n(),
mean_Slope = mean(Slope, na.rm = TRUE),
.groups = "drop"
)
cat("Number of Blocks after averaging:", nrow(block_avg_slope), "\n")
## Number of Blocks after averaging: 264
head(block_avg_slope)
## # A tibble: 6 × 8
## Block District Agzone Officer cell Quali n_seasons mean_Slope
## <chr> <chr> <chr> <chr> <chr> <fct> <int> <dbl>
## 1 102_10 Bugarama Bugarama Wellars Hangabashi Very Low 6 15
## 2 102_11 Bugarama Bugarama Wellars Hangabashi Very Low 6 15
## 3 102_12 Bugarama Bugarama Wellars Hangabashi Very High 6 15
## 4 102_13 Bugarama Bugarama Wellars Hangabashi Very Low 6 20
## 5 102_14 Bugarama Bugarama Wellars Hangabashi Very High 6 12
## 6 102_2 Bugarama Bugarama Wellars Hangabashi Very High 6 14
# --- Step 2: Tukey HSD across Quali groups ---
model <- aov(mean_Slope ~ Quali, data = block_avg_slope)
tukey <- TukeyHSD(model)
letters_df <- multcompLetters4(model, tukey)$Quali$Letters
letters_tbl <- data.frame(Quali = names(letters_df), letter = letters_df, row.names = NULL)
# --- Step 3: N, mean, y-position per Quali group ---
label_stats <- block_avg_slope %>%
group_by(Quali) %>%
summarise(
n = n(),
mean_v = mean(mean_Slope, na.rm = TRUE),
y_max = max(mean_Slope, na.rm = TRUE),
.groups = "drop"
)
label_data <- label_stats %>%
left_join(letters_tbl, by = "Quali") %>%
mutate(
label = paste0(letter, "\nn=", n, "\nmean=", round(mean_v, 2)),
y_pos = y_max * 1.05
)
# --- Step 4: Plot ---
p_block_avg_slope <- ggplot(block_avg_slope, aes(x = Quali, y = mean_Slope, fill = Quali)) +
geom_boxplot(na.rm = TRUE, outlier.shape = NA, alpha = 0.6) +
geom_jitter(width = 0.15, height = 0, alpha = 0.6, size = 1.8, color = "black") +
geom_text(data = label_data,
aes(x = Quali, y = y_pos, label = label),
inherit.aes = FALSE, vjust = 0, size = 3, fontface = "bold", lineheight = 0.85) +
labs(title = "Average Slope per Block (across seasons) by Quali",
x = "Quality of the trial", y = "Mean Slope (per Block)") +
theme_minimal() +
theme(axis.text.x = element_text(angle = 45, hjust = 1), legend.position = "none") +
scale_y_continuous(expand = expansion(mult = c(0.05, 0.2)))
print(p_block_avg_slope)
ggsave(file.path(output_dir, "boxplot_Slope_by_Block_avg_Quali_dots.png"), p_block_avg_slope, width = 8, height = 6)
Effect on position on th eslope?
library(dplyr)
library(ggplot2)
library(tidyr)
# --- Step 0: add IQF_position_slope into check_CA1 ---
yield_subset <- yield_CA1 %>%
select(Block, IQR_Season, IQF_environment, IQR_sex_HH_head,
ICF_Lat, ICF_Long, ICF_Alt, Y_ratio, env_index, yield_A, yield_C,
Slope, IQF_position_slope)
check_CA1 <- sel_subset %>%
left_join(yield_subset, by = "Block") %>%
mutate(
data_flag = ifelse(is.na(IQR_Season), "no_yield_data_found", "ok"),
# clean junk codes into real NA
IQF_position_slope = na_if(IQF_position_slope, "---"),
IQF_position_slope = na_if(IQF_position_slope, "#N/A")
)
## Warning in left_join(., yield_subset, by = "Block"): Detected an unexpected many-to-many relationship between `x` and `y`.
## ℹ Row 1 of `x` matches multiple rows in `y`.
## ℹ Row 1855 of `y` matches multiple rows in `x`.
## ℹ If a many-to-many relationship is expected, set `relationship =
## "many-to-many"` to silence this warning.
check_CA1_clean <- check_CA1 %>%
filter(data_flag != "no_yield_data_found") %>%
mutate(Quali = factor(Quali, levels = c("Very Low", "Low", "Medium", "High", "Very High")))
# --- Step 1: most common slope-position per Block across seasons ---
get_mode <- function(x) {
x <- x[!is.na(x)]
if (length(x) == 0) return(NA_character_)
names(sort(table(x), decreasing = TRUE))[1]
}
block_slope_pos <- check_CA1_clean %>%
group_by(Block, District, Agzone, Officer, cell, Quali) %>%
summarise(
n_seasons = n(),
slope_pos_mode = get_mode(IQF_position_slope),
.groups = "drop"
) %>%
filter(!is.na(slope_pos_mode))
cat("Number of Blocks with slope-position data:", nrow(block_slope_pos), "\n")
## Number of Blocks with slope-position data: 261
# --- Step 2: Chi-square test of independence (Quali x slope position) ---
cont_table <- table(block_slope_pos$Quali, block_slope_pos$slope_pos_mode)
print(cont_table)
##
## Backslope Flat Footslope Shoulder Summit
## Very Low 45 3 9 2 7
## Low 6 0 0 3 1
## Medium 26 3 2 5 3
## High 4 0 0 1 0
## Very High 93 11 18 16 3
chisq_result <- chisq.test(cont_table)
## Warning in chisq.test(cont_table): Chi-squared approximation may be incorrect
print(chisq_result)
##
## Pearson's Chi-squared test
##
## data: cont_table
## X-squared = 20.858, df = 16, p-value = 0.184
# --- Step 3: proportions for plotting ---
prop_data <- block_slope_pos %>%
count(Quali, slope_pos_mode) %>%
group_by(Quali) %>%
mutate(prop = n / sum(n)) %>%
ungroup()
# --- Step 4: stacked bar chart ---
p_slope_pos <- ggplot(prop_data, aes(x = Quali, y = prop, fill = slope_pos_mode)) +
geom_col(position = "stack", color = "white") +
geom_text(aes(label = ifelse(prop > 0.04, paste0(round(prop*100), "%"), "")),
position = position_stack(vjust = 0.5), size = 3, color = "black") +
labs(title = paste0("Slope position by Quali (Chi-sq p = ",
signif(chisq_result$p.value, 3), ")"),
x = "Quality of the trial", y = "Proportion of Blocks", fill = "Slope position") +
theme_minimal() +
theme(axis.text.x = element_text(angle = 45, hjust = 1))
print(p_slope_pos)
ggsave(file.path(output_dir, "barplot_slopeposition_by_Quali.png"), p_slope_pos, width = 8, height = 6)
Map of plot sdistribution
library(dplyr)
library(ggplot2)
# --- Step 1: one Lat/Long per Block ---
block_coords <- check_CA1_clean %>%
group_by(Block, Quali) %>%
summarise(
Lat = mean(ICF_Lat, na.rm = TRUE),
Long = mean(ICF_Long, na.rm = TRUE),
.groups = "drop"
) %>%
filter(!is.na(Lat), !is.na(Long))
cat("Number of Blocks with coordinates:", nrow(block_coords), "\n")
## Number of Blocks with coordinates: 264
# --- Step 2: get Rwanda country boundary ---
# install.packages("rnaturalearth") # run once if not installed
# install.packages("sf") # run once if not installed
# install.packages("geodata") # run once if not installed
library(geodata)
## 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
library(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 TRUE
library(dplyr)
library(ggplot2)
# --- Step 1: Rwanda district boundaries ---
rwanda_district <- gadm(country = "RWA", level = 2, path = tempdir()) %>% st_as_sf()
# --- Step 2: Block coordinates ---
block_coords <- check_CA1_clean %>%
group_by(Block, Quali) %>%
summarise(
Lat = mean(ICF_Lat, na.rm = TRUE),
Long = mean(ICF_Long, na.rm = TRUE),
.groups = "drop"
) %>%
filter(!is.na(Lat), !is.na(Long)) %>%
mutate(Quali = factor(Quali, levels = c("Very Low", "Low", "Medium", "High", "Very High")))
# --- Step 3: Plot with district boundaries + smaller, jittered points ---
set.seed(123) # for reproducible jitter placement
p_map <- ggplot() +
geom_sf(data = rwanda_district, fill = "gray97", color = "black", linewidth = 0.5) +
geom_jitter(data = block_coords, aes(x = Long, y = Lat, color = Quali),
width = 0.02, height = 0.02, size = 1.2, alpha = 0.8) +
scale_color_manual(values = c("Very Low" = "#d73027", "Low" = "#fc8d59",
"Medium" = "#fee08b", "High" = "#91cf60",
"Very High" = "#1a9850")) +
coord_sf(xlim = c(28.8, 30.9), ylim = c(-2.85, -1.0), expand = FALSE) +
labs(title = "Farm Blocks by Quali, with District boundaries",
x = NULL, y = NULL, color = "Quali") +
theme_minimal() +
theme(panel.grid = element_blank())
print(p_map)
ggsave(file.path(output_dir, "map_Blocks_Quali_district_boundaries.png"), p_map, width = 10, height = 10, dpi = 300)