Clean_Results
Diego
2026-02-18
# Load packages
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)
CONV_CA_Ratio <- read_csv(
"D:/Mes Donnees/OAF_CIRAD/CA_MoU/data/processed/CONV_CA_Ratio.csv",
col_types = cols(
IQR_weed1_complexity = col_character()
)
)
drivers_data <- read_csv("data/processed/drivers_data.csv")
r1 <- read_csv("data/processed/r1.csv")
library(dplyr)Mix model to evaluate effect of crop, season, year and interactions on CA:CONV ratio
Description of the analysis: Linear mixed-effects model fitted to the log-transformed yield ratio (log(Y_ratio + 1)) comparing Conservation Agriculture (CA) to Conventional Agriculture (CONV). The model included:Fixed effects: crop, year, IQF_environment, and all two- and three-way interactions Random effect: IQR_block The treatment variable (Treat_code) was excluded in this model (fit2_conv) Model performance was compared to a simpler model (fit1_conv) using BIC. The summary of fixed effects and significance tests were obtained via ANOVA with Satterthwaite approximation.
Description of the results: On the analysis of the log-transformed CA:CONV yield ratio, there was a significant effect of crop, year, and environment, as well as their two- and three-way interactions.
# Model 1: no log transform, no Treat_code
fit1_conv <- lmer(
Y_ratio ~ crop + year + IQR_Season + IQF_agzone +
IQR_Season:IQF_agzone +
(1 | IQR_block),
data = CONV_CA_Ratio
)
# Model 2: log transform, no Treat_code
fit2_conv <- lmer(
log(Y_ratio + 1) ~ crop + year + crop:year + IQF_agzone +
crop:IQF_agzone +
year:IQF_agzone +
year:crop:IQF_agzone +
year:crop +
(1 | IQR_block),
data = CONV_CA_Ratio
)
# Model comparison
BIC(fit1_conv, fit2_conv)# View model summaries
summary(fit1_conv)## Linear mixed model fit by REML ['lmerMod']
## Formula:
## Y_ratio ~ crop + year + IQR_Season + IQF_agzone + IQR_Season:IQF_agzone +
## (1 | IQR_block)
## Data: CONV_CA_Ratio
##
## REML criterion at convergence: 1851.9
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -2.6225 -0.5398 -0.0952 0.3766 9.8431
##
## Random effects:
## Groups Name Variance Std.Dev.
## IQR_block (Intercept) 0.01578 0.1256
## Residual 0.09353 0.3058
## Number of obs: 2818, groups: IQR_block, 565
##
## Fixed effects:
## Estimate Std. Error t value
## (Intercept) 6.2252926 0.8046305 7.737
## cropMaize -0.1538182 0.0685580 -2.244
## year -0.2198388 0.0321626 -6.835
## IQR_Season23B -0.2252049 0.0901646 -2.498
## IQR_Season24A -0.2275483 0.0540685 -4.209
## IQR_Season24B -0.3812046 0.0722806 -5.274
## IQF_agzoneBugarama -0.0693830 0.0744988 -0.931
## IQF_agzoneCentral Plateau -0.0635196 0.0578590 -1.098
## IQF_agzoneCongo Nile -0.1497084 0.0605970 -2.471
## IQF_agzoneCyangugu -0.0660860 0.0618429 -1.069
## IQF_agzoneEastern Ridges -0.0871046 0.0601740 -1.448
## IQF_agzoneEastern Savanah -0.2385440 0.0635084 -3.756
## IQF_agzoneLake Kivu -0.1564174 0.0608275 -2.571
## IQF_agzoneMayaga-Bugesera 0.0039428 0.0626944 0.063
## IQF_agzoneVolcanic cones -0.2072991 0.0778891 -2.661
## IQR_Season23B:IQF_agzoneBugarama -0.0597157 0.1018381 -0.586
## IQR_Season24A:IQF_agzoneBugarama 0.4370113 0.1032515 4.232
## IQR_Season24B:IQF_agzoneBugarama 0.0879833 0.1034762 0.850
## IQR_Season25A:IQF_agzoneBugarama 0.4689238 0.1080181 4.341
## IQR_Season25B:IQF_agzoneBugarama 0.2178256 0.1064512 2.046
## IQR_Season23B:IQF_agzoneCentral Plateau -0.0123275 0.0771825 -0.160
## IQR_Season24A:IQF_agzoneCentral Plateau 0.3660999 0.0784522 4.667
## IQR_Season24B:IQF_agzoneCentral Plateau 0.3197537 0.0788915 4.053
## IQR_Season25A:IQF_agzoneCentral Plateau 0.3156120 0.0827658 3.813
## IQR_Season25B:IQF_agzoneCentral Plateau 0.3812980 0.0817404 4.665
## IQR_Season23B:IQF_agzoneCongo Nile 0.0897171 0.0830668 1.080
## IQR_Season24A:IQF_agzoneCongo Nile 0.1987473 0.0840807 2.364
## IQR_Season24B:IQF_agzoneCongo Nile 0.1047993 0.0836076 1.253
## IQR_Season25A:IQF_agzoneCongo Nile 0.1646176 0.0870750 1.891
## IQR_Season25B:IQF_agzoneCongo Nile 0.3285289 0.0865113 3.798
## IQR_Season23B:IQF_agzoneCyangugu 0.0265940 0.0824284 0.323
## IQR_Season24A:IQF_agzoneCyangugu 0.2248523 0.0879746 2.556
## IQR_Season24B:IQF_agzoneCyangugu 0.2169128 0.0852752 2.544
## IQR_Season25A:IQF_agzoneCyangugu 0.2225999 0.0896337 2.483
## IQR_Season25B:IQF_agzoneCyangugu 0.2116234 0.0882016 2.399
## IQR_Season23B:IQF_agzoneEastern Ridges -0.2841609 0.0809814 -3.509
## IQR_Season24A:IQF_agzoneEastern Ridges 0.3524881 0.0820912 4.294
## IQR_Season24B:IQF_agzoneEastern Ridges 0.3256443 0.0825896 3.943
## IQR_Season25A:IQF_agzoneEastern Ridges 0.2472866 0.0862849 2.866
## IQR_Season25B:IQF_agzoneEastern Ridges 0.1001398 0.0850444 1.177
## IQR_Season23B:IQF_agzoneEastern Savanah -0.0647137 0.0841196 -0.769
## IQR_Season24A:IQF_agzoneEastern Savanah 0.4032337 0.0883363 4.565
## IQR_Season24B:IQF_agzoneEastern Savanah 0.2706519 0.0878230 3.082
## IQR_Season25A:IQF_agzoneEastern Savanah 0.4123700 0.0910739 4.528
## IQR_Season25B:IQF_agzoneEastern Savanah 0.4184354 0.0908518 4.606
## IQR_Season23B:IQF_agzoneLake Kivu 0.0005337 0.0820111 0.007
## IQR_Season24A:IQF_agzoneLake Kivu 0.4055660 0.0836467 4.849
## IQR_Season24B:IQF_agzoneLake Kivu 0.3738662 0.0839369 4.454
## IQR_Season25A:IQF_agzoneLake Kivu 0.3787711 0.0891277 4.250
## IQR_Season25B:IQF_agzoneLake Kivu 0.2912298 0.0868787 3.352
## IQR_Season23B:IQF_agzoneMayaga-Bugesera -0.1557215 0.0840599 -1.853
## IQR_Season24A:IQF_agzoneMayaga-Bugesera 0.4127703 0.0842678 4.898
## IQR_Season24B:IQF_agzoneMayaga-Bugesera 0.3302302 0.0850138 3.884
## IQR_Season25A:IQF_agzoneMayaga-Bugesera 0.3571996 0.0898285 3.976
## IQR_Season25B:IQF_agzoneMayaga-Bugesera 0.1919620 0.0878959 2.184
## IQR_Season23B:IQF_agzoneVolcanic cones 0.1060746 0.1017404 1.043
## IQR_Season24A:IQF_agzoneVolcanic cones 0.3292551 0.1075476 3.061
## IQR_Season24B:IQF_agzoneVolcanic cones 0.4269823 0.1038557 4.111
## IQR_Season25A:IQF_agzoneVolcanic cones 0.4333492 0.1067415 4.060
## IQR_Season25B:IQF_agzoneVolcanic cones 0.4357012 0.1065441 4.089
## fit warnings:
## fixed-effect model matrix is rank deficient so dropping 2 columns / coefficientssummary(fit2_conv)## Linear mixed model fit by REML ['lmerMod']
## Formula:
## log(Y_ratio + 1) ~ crop + year + crop:year + IQF_agzone + crop:IQF_agzone +
## year:IQF_agzone + year:crop:IQF_agzone + year:crop + (1 | IQR_block)
## Data: CONV_CA_Ratio
##
## REML criterion at convergence: -1919.1
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -3.0632 -0.5634 -0.0533 0.5013 6.6506
##
## Random effects:
## Groups Name Variance Std.Dev.
## IQR_block (Intercept) 0.004464 0.06681
## Residual 0.024142 0.15538
## Number of obs: 2818, groups: IQR_block, 565
##
## Fixed effects:
## Estimate Std. Error t value
## (Intercept) 2.127314 0.388137 5.481
## cropMaize 1.569377 0.545271 2.878
## year -0.066057 0.016223 -4.072
## IQF_agzoneBugarama -1.863768 0.666638 -2.796
## IQF_agzoneCentral Plateau -2.570038 0.501295 -5.127
## IQF_agzoneCongo Nile -1.456010 0.543418 -2.679
## IQF_agzoneCyangugu -1.081076 0.540270 -2.001
## IQF_agzoneEastern Ridges -2.593732 0.525743 -4.933
## IQF_agzoneEastern Savanah -3.356843 0.548409 -6.121
## IQF_agzoneLake Kivu -1.846136 0.537726 -3.433
## IQF_agzoneMayaga-Bugesera -2.279403 0.541530 -4.209
## IQF_agzoneVolcanic cones -2.293946 0.636451 -3.604
## cropMaize:year -0.066244 0.022804 -2.905
## cropMaize:IQF_agzoneBugarama -1.371371 0.928541 -1.477
## cropMaize:IQF_agzoneCentral Plateau 0.345749 0.705232 0.490
## cropMaize:IQF_agzoneCongo Nile 0.251533 0.753019 0.334
## cropMaize:IQF_agzoneCyangugu -0.677607 0.761520 -0.890
## cropMaize:IQF_agzoneEastern Ridges 0.706534 0.737805 0.958
## cropMaize:IQF_agzoneEastern Savanah 0.406712 0.773141 0.526
## cropMaize:IQF_agzoneLake Kivu -0.776177 0.756483 -1.026
## cropMaize:IQF_agzoneMayaga-Bugesera -0.083462 0.764416 -0.109
## cropMaize:IQF_agzoneVolcanic cones -0.872411 0.905554 -0.963
## year:IQF_agzoneBugarama 0.078109 0.027812 2.808
## year:IQF_agzoneCentral Plateau 0.110758 0.020933 5.291
## year:IQF_agzoneCongo Nile 0.060600 0.022672 2.673
## year:IQF_agzoneCyangugu 0.046836 0.022573 2.075
## year:IQF_agzoneEastern Ridges 0.106688 0.021944 4.862
## year:IQF_agzoneEastern Savanah 0.138989 0.022929 6.062
## year:IQF_agzoneLake Kivu 0.078297 0.022454 3.487
## year:IQF_agzoneMayaga-Bugesera 0.097948 0.022598 4.334
## year:IQF_agzoneVolcanic cones 0.098347 0.026564 3.702
## cropMaize:year:IQF_agzoneBugarama 0.062517 0.038790 1.612
## cropMaize:year:IQF_agzoneCentral Plateau -0.013838 0.029464 -0.470
## cropMaize:year:IQF_agzoneCongo Nile -0.010720 0.031448 -0.341
## cropMaize:year:IQF_agzoneCyangugu 0.028755 0.031836 0.903
## cropMaize:year:IQF_agzoneEastern Ridges -0.024853 0.030820 -0.806
## cropMaize:year:IQF_agzoneEastern Savanah -0.014654 0.032313 -0.453
## cropMaize:year:IQF_agzoneLake Kivu 0.033737 0.031622 1.067
## cropMaize:year:IQF_agzoneMayaga-Bugesera 0.006606 0.031928 0.207
## cropMaize:year:IQF_agzoneVolcanic cones 0.035376 0.037777 0.936# Fixed effects significance
anova(fit2_conv)# Residual diagnostics
plot(fit2_conv)
qqnorm(resid(fit2_conv)); qqline(resid(fit2_conv))
library(lmerTest) # masks lme4::lmer automatically
# Re-run your models exactly as before — p-values now appear in summary()
fit2_conv <- lmer(
log(Y_ratio + 1) ~ crop + year + crop:year + IQF_agzone +
crop:IQF_agzone + year:IQF_agzone + year:crop:IQF_agzone +
(1 | IQR_block),
data = CONV_CA_Ratio
)
summary(fit2_conv) # now shows p-values in fixed effects table## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula:
## log(Y_ratio + 1) ~ crop + year + crop:year + IQF_agzone + crop:IQF_agzone +
## year:IQF_agzone + year:crop:IQF_agzone + (1 | IQR_block)
## Data: CONV_CA_Ratio
##
## REML criterion at convergence: -1919.1
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -3.0632 -0.5634 -0.0533 0.5013 6.6506
##
## Random effects:
## Groups Name Variance Std.Dev.
## IQR_block (Intercept) 0.004464 0.06681
## Residual 0.024142 0.15538
## Number of obs: 2818, groups: IQR_block, 565
##
## Fixed effects:
## Estimate Std. Error df
## (Intercept) 2.127e+00 3.881e-01 2.330e+03
## cropMaize 1.569e+00 5.453e-01 2.292e+03
## year -6.606e-02 1.622e-02 2.330e+03
## IQF_agzoneBugarama -1.864e+00 6.666e-01 2.333e+03
## IQF_agzoneCentral Plateau -2.570e+00 5.013e-01 2.315e+03
## IQF_agzoneCongo Nile -1.456e+00 5.434e-01 2.343e+03
## IQF_agzoneCyangugu -1.081e+00 5.403e-01 2.332e+03
## IQF_agzoneEastern Ridges -2.594e+00 5.257e-01 2.308e+03
## IQF_agzoneEastern Savanah -3.357e+00 5.484e-01 2.351e+03
## IQF_agzoneLake Kivu -1.846e+00 5.377e-01 2.338e+03
## IQF_agzoneMayaga-Bugesera -2.279e+00 5.415e-01 2.307e+03
## IQF_agzoneVolcanic cones -2.294e+00 6.365e-01 2.306e+03
## cropMaize:year -6.624e-02 2.280e-02 2.290e+03
## cropMaize:IQF_agzoneBugarama -1.371e+00 9.285e-01 2.309e+03
## cropMaize:IQF_agzoneCentral Plateau 3.457e-01 7.052e-01 2.293e+03
## cropMaize:IQF_agzoneCongo Nile 2.515e-01 7.530e-01 2.307e+03
## cropMaize:IQF_agzoneCyangugu -6.776e-01 7.615e-01 2.305e+03
## cropMaize:IQF_agzoneEastern Ridges 7.065e-01 7.378e-01 2.295e+03
## cropMaize:IQF_agzoneEastern Savanah 4.067e-01 7.731e-01 2.326e+03
## cropMaize:IQF_agzoneLake Kivu -7.762e-01 7.565e-01 2.306e+03
## cropMaize:IQF_agzoneMayaga-Bugesera -8.346e-02 7.644e-01 2.302e+03
## cropMaize:IQF_agzoneVolcanic cones -8.724e-01 9.056e-01 2.315e+03
## year:IQF_agzoneBugarama 7.811e-02 2.781e-02 2.329e+03
## year:IQF_agzoneCentral Plateau 1.108e-01 2.093e-02 2.314e+03
## year:IQF_agzoneCongo Nile 6.060e-02 2.267e-02 2.343e+03
## year:IQF_agzoneCyangugu 4.684e-02 2.257e-02 2.332e+03
## year:IQF_agzoneEastern Ridges 1.067e-01 2.194e-02 2.306e+03
## year:IQF_agzoneEastern Savanah 1.390e-01 2.293e-02 2.353e+03
## year:IQF_agzoneLake Kivu 7.830e-02 2.245e-02 2.337e+03
## year:IQF_agzoneMayaga-Bugesera 9.795e-02 2.260e-02 2.305e+03
## year:IQF_agzoneVolcanic cones 9.835e-02 2.656e-02 2.304e+03
## cropMaize:year:IQF_agzoneBugarama 6.252e-02 3.879e-02 2.307e+03
## cropMaize:year:IQF_agzoneCentral Plateau -1.384e-02 2.946e-02 2.292e+03
## cropMaize:year:IQF_agzoneCongo Nile -1.072e-02 3.145e-02 2.306e+03
## cropMaize:year:IQF_agzoneCyangugu 2.875e-02 3.184e-02 2.304e+03
## cropMaize:year:IQF_agzoneEastern Ridges -2.485e-02 3.082e-02 2.294e+03
## cropMaize:year:IQF_agzoneEastern Savanah -1.465e-02 3.231e-02 2.325e+03
## cropMaize:year:IQF_agzoneLake Kivu 3.374e-02 3.162e-02 2.305e+03
## cropMaize:year:IQF_agzoneMayaga-Bugesera 6.606e-03 3.193e-02 2.301e+03
## cropMaize:year:IQF_agzoneVolcanic cones 3.538e-02 3.778e-02 2.314e+03
## t value Pr(>|t|)
## (Intercept) 5.481 4.69e-08 ***
## cropMaize 2.878 0.004037 **
## year -4.072 4.82e-05 ***
## IQF_agzoneBugarama -2.796 0.005220 **
## IQF_agzoneCentral Plateau -5.127 3.19e-07 ***
## IQF_agzoneCongo Nile -2.679 0.007428 **
## IQF_agzoneCyangugu -2.001 0.045509 *
## IQF_agzoneEastern Ridges -4.933 8.65e-07 ***
## IQF_agzoneEastern Savanah -6.121 1.09e-09 ***
## IQF_agzoneLake Kivu -3.433 0.000607 ***
## IQF_agzoneMayaga-Bugesera -4.209 2.66e-05 ***
## IQF_agzoneVolcanic cones -3.604 0.000320 ***
## cropMaize:year -2.905 0.003708 **
## cropMaize:IQF_agzoneBugarama -1.477 0.139836
## cropMaize:IQF_agzoneCentral Plateau 0.490 0.623995
## cropMaize:IQF_agzoneCongo Nile 0.334 0.738386
## cropMaize:IQF_agzoneCyangugu -0.890 0.373661
## cropMaize:IQF_agzoneEastern Ridges 0.958 0.338357
## cropMaize:IQF_agzoneEastern Savanah 0.526 0.598902
## cropMaize:IQF_agzoneLake Kivu -1.026 0.304983
## cropMaize:IQF_agzoneMayaga-Bugesera -0.109 0.913067
## cropMaize:IQF_agzoneVolcanic cones -0.963 0.335447
## year:IQF_agzoneBugarama 2.808 0.005020 **
## year:IQF_agzoneCentral Plateau 5.291 1.33e-07 ***
## year:IQF_agzoneCongo Nile 2.673 0.007572 **
## year:IQF_agzoneCyangugu 2.075 0.038113 *
## year:IQF_agzoneEastern Ridges 4.862 1.24e-06 ***
## year:IQF_agzoneEastern Savanah 6.062 1.56e-09 ***
## year:IQF_agzoneLake Kivu 3.487 0.000497 ***
## year:IQF_agzoneMayaga-Bugesera 4.334 1.52e-05 ***
## year:IQF_agzoneVolcanic cones 3.702 0.000219 ***
## cropMaize:year:IQF_agzoneBugarama 1.612 0.107172
## cropMaize:year:IQF_agzoneCentral Plateau -0.470 0.638636
## cropMaize:year:IQF_agzoneCongo Nile -0.341 0.733227
## cropMaize:year:IQF_agzoneCyangugu 0.903 0.366503
## cropMaize:year:IQF_agzoneEastern Ridges -0.806 0.420107
## cropMaize:year:IQF_agzoneEastern Savanah -0.453 0.650239
## cropMaize:year:IQF_agzoneLake Kivu 1.067 0.286141
## cropMaize:year:IQF_agzoneMayaga-Bugesera 0.207 0.836113
## cropMaize:year:IQF_agzoneVolcanic cones 0.936 0.349149
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1anova(fit2_conv) # now shows F-stats + p-values with dfDescription of the performance of CA across Environments and Districts, and from season 1 to 3
Description of the analysis: Used eemeans from the mixed-effects model to summarize the overall CA:CONV yield ratio, and to assess the effects of crop, year, and environment. Results were back-transformed for interpretation on the original yield ratio scale. Interaction effects were explored and visualized.
Description of the results: Conventional outperforms CA as a whole experimental average; overall CA:CONV ratio = 0.787 (eemeans for the same model shown above). Advantage of CONV was similar for maize and beans (t-tests of individual parameter estimates p=0.7). Although there was a significant effect of year, there was no trend for CA to perform better from year 1 to 3. Also, there was a clear effect of the environment with transition and mid-low-wet being more favourable for CA
Overal difference
overall_emm <- emmeans(fit2_conv, ~ 1, type = "response")
print(overall_emm)## 1 response SE df lower.CL upper.CL
## overall 0.783 0.00772 528 0.768 0.798
##
## Results are averaged over the levels of: crop, IQF_agzone
## Degrees-of-freedom method: kenward-roger
## Confidence level used: 0.95
## Intervals are back-transformed from the log(mu + 1) scaleeemeans to interpret the effects of categorical predictors and their interactions
library(emmeans)
# Compute emmeans for crop within year × environment
emm_crop_year_env <- emmeans(fit2_conv, ~ crop | year * IQF_agzone, type = "response")
emm_df <- as.data.frame(emm_crop_year_env)eemeans per crop
emm_df %>%
group_by(crop) %>%
summarise(avg_ratio = mean(response, na.rm = TRUE))Does the % wins of CA increse progresibely from year 1 to 3? Separated by crop
library(dplyr)
library(ggplot2)
library(broom)
library(tidyr)
# 1. Prepare data
trend_data <- CONV_CA_Ratio %>%
filter(
is.finite(Y_ratio),
year %in% c("23", "24", "25"),
IQR_SeasAB %in% c("A", "B")
) %>%
mutate(
year = as.numeric(year),
crop_type = ifelse(IQR_SeasAB == "A", "Maize", "Beans")
)
# 2. Function to compute slope counts
get_slope_label <- function(df) {
# One slope per block
slopes <- df %>%
group_by(IQR_block) %>%
filter(n_distinct(year) >= 2) %>%
summarise(slope = coef(lm(Y_ratio ~ year))[2], .groups = "drop") %>%
mutate(
direction = case_when(
slope > 0 ~ "↑ Positive",
slope < 0 ~ "↓ Negative",
TRUE ~ "= Flat"
)
) %>%
count(direction) %>%
complete(direction, fill = list(n = 0))
# Format label text
paste0(
"↑ Positive: ", slopes$n[slopes$direction == "↑ Positive"], "\n",
"↓ Negative: ", slopes$n[slopes$direction == "↓ Negative"], "\n",
"= Flat: ", slopes$n[slopes$direction == "= Flat"]
)
}
# Labels
label_maize <- get_slope_label(trend_data %>% filter(crop_type == "Maize"))
label_beans <- get_slope_label(trend_data %>% filter(crop_type == "Beans"))
# Common y-limit helper
y_max <- max(trend_data$Y_ratio, na.rm = TRUE)
# 3. MAIZE plot
p_maize <- ggplot(
trend_data %>% filter(crop_type == "Maize"),
aes(x = year, y = Y_ratio)
) +
geom_point(alpha = 0.35, size = 2) +
# Block-level linear trends
geom_smooth(
aes(group = IQR_block),
method = "lm",
se = FALSE,
color = "steelblue",
linewidth = 0.6
) +
# Global trend
geom_smooth(
method = "lm",
se = FALSE,
color = "red",
linewidth = 1.3
) +
# Annotation
annotate(
"label",
x = 23.05,
y = y_max * 0.95,
label = label_maize,
hjust = 0,
vjust = 1,
size = 4.5,
fontface = "bold",
fill = "white",
label.size = 0.4
) +
scale_x_continuous(breaks = c(23, 24, 25)) +
labs(
title = "Maize",
x = "Year",
y = "Yield Ratio (CA:CT)"
) +
theme_minimal() +
theme(
plot.title = element_text(face = "bold", hjust = 0.5, size = 16)
)
# 4. BEANS plot
p_beans <- ggplot(
trend_data %>% filter(crop_type == "Beans"),
aes(x = year, y = Y_ratio)
) +
geom_point(alpha = 0.35, size = 2) +
# Block-level linear trends
geom_smooth(
aes(group = IQR_block),
method = "lm",
se = FALSE,
color = "darkorange",
linewidth = 0.6
) +
# Global trend
geom_smooth(
method = "lm",
se = FALSE,
color = "red",
linewidth = 1.3
) +
# Annotation
annotate(
"label",
x = 23.05,
y = y_max * 0.95,
label = label_beans,
hjust = 0,
vjust = 1,
size = 4.5,
fontface = "bold",
fill = "white",
label.size = 0.4
) +
scale_x_continuous(breaks = c(23, 24, 25)) +
labs(
title = "Beans",
x = "Year",
y = "Yield Ratio (CA:CT)"
) +
theme_minimal() +
theme(
plot.title = element_text(face = "bold", hjust = 0.5, size = 16)
)
# 5. Print plots
print(p_maize)
print(p_beans)
Lets define what is an acceptable yiel dloss
Based on the relation among farmers “preference” and relative yield of the season s
#Figure 4A
library(dplyr)
library(minpack.lm)
library(nlraa)
# --- Sort and group ────────────────────────────────────────────────────────────
CONV_CA_Ratio <- CONV_CA_Ratio[order(CONV_CA_Ratio$Y_ratio), ]
CONV_CA_Ratio$group15 <- rep(1:ceiling(nrow(CONV_CA_Ratio)/20), each = 20)[1:nrow(CONV_CA_Ratio)]
grouped_data <- CONV_CA_Ratio %>%
group_by(group15) %>%
summarise(
mean_Y_ratio = mean(Y_ratio, na.rm = TRUE),
mean_DQ = mean(DQ_future_adoption, na.rm = TRUE),
.groups = "drop"
) %>%
filter(mean_Y_ratio <= 2.5) # ← remove groups with mean > 2.5
# --- Fit linear-plateau model ─────────────────────────────────────────────────
fit_lp <- nlsLM(
mean_DQ ~ SSlinp(mean_Y_ratio, a, b, xs),
data = grouped_data,
start = list(a = 0.3, b = 0.5, xs = 1.0),
lower = c(a = 0, b = 0, xs = min(grouped_data$mean_Y_ratio)),
upper = c(a = 1, b = 5, xs = max(grouped_data$mean_Y_ratio))
)
summary(fit_lp)##
## Formula: mean_DQ ~ SSlinp(mean_Y_ratio, a, b, xs)
##
## Parameters:
## Estimate Std. Error t value Pr(>|t|)
## a 0.16995 0.03443 4.936 2.28e-06 ***
## b 0.58563 0.04568 12.822 < 2e-16 ***
## xs 1.13279 0.05886 19.244 < 2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.1116 on 137 degrees of freedom
##
## Number of iterations to convergence: 10
## Achieved convergence tolerance: 1.49e-08# --- Extract parameters ────────────────────────────────────────────────────────
a_est <- coef(fit_lp)["a"]
b_est <- coef(fit_lp)["b"]
xs_est <- coef(fit_lp)["xs"]
plateau <- a_est + b_est * xs_est
# --- R² ───────────────────────────────────────────────────────────────────────
ss_tot <- sum((grouped_data$mean_DQ - mean(grouped_data$mean_DQ))^2)
R2 <- 1 - sum(residuals(fit_lp)^2) / ss_tot
cat("Breakpoint =", round(xs_est, 3), "\n")## Breakpoint = 1.133cat("Plateau =", round(plateau, 3), "\n")## Plateau = 0.833cat("R² =", round(R2, 3), "\n")## R² = 0.632# --- Target: 15% below plateau ────────────────────────────────────────────────
target <- plateau - 0.15
Y_ratio_target <- (target - a_est) / b_est
cat("Y_ratio for target =", round(Y_ratio_target, 3), "\n")## Y_ratio for target = 0.877# --- Prediction grid ──────────────────────────────────────────────────────────
newdata <- data.frame(
mean_Y_ratio = seq(min(grouped_data$mean_Y_ratio),
max(grouped_data$mean_Y_ratio),
length.out = 500)
)
newdata$pred <- predict(fit_lp, newdata)
# --- Plot ─────────────────────────────────────────────────────────────────────
plot(grouped_data$mean_Y_ratio,
grouped_data$mean_DQ,
pch = 16,
xlab = "Mean CA:CT Yield Ratio (groups of 20)",
ylab = "Proportion of farmers reporting future adoption",
ylim = c(0, 1))
# --- Panel label ──────────────────────────────────────────────────────────────
#mtext("A", side = 3, adj = 0, line = 0.5, cex = 1.5, font = 2) # ← top left
# --- Panel label inside ───────────────────────────────────────────────────────
text(x = par("usr")[1], # ← left edge of plot area
y = par("usr")[4], # ← top edge of plot area
labels = "A",
adj = c(-0.3, 1.3), # ← small offset from corner
cex = 1.8, font = 2) # ← size and bold
lines(newdata$mean_Y_ratio, newdata$pred,
col = "blue", lwd = 2)
# --- Breakpoint lines ─────────────────────────────────────────────────────────
abline(v = xs_est, col = "grey", lwd = 3, lty = 2)
abline(h = plateau, col = "grey", lwd = 2, lty = 3)
points(xs_est, plateau, pch = 19, col = "grey")
# --- Target lines ─────────────────────────────────────────────────────────────
abline(h = target, col = "#009900", lwd = 3, lty = 2)
abline(v = Y_ratio_target, col = "#009900", lwd = 3, lty = 1)
points(Y_ratio_target, target, pch = 19, col = "#009900", cex = 1.5)
# --- Annotations ──────────────────────────────────────────────────────────────
text(xs_est, 0.0,
labels = paste0("Breakpoint = ", round(xs_est, 3)),
col = "black", pos = 4)
text(xs_est, plateau + 0.05,
labels = paste0("Adoption at Plateau = ", round(plateau, 3)),
col = "black", pos = 4)
text(Y_ratio_target, 0.2,
labels = paste0("Yield Ratio (Plateau − 0.15) = ", round(Y_ratio_target, 3)),
col = "#009900", pos = 4)
text(max(grouped_data$mean_Y_ratio), 0.95,
labels = paste0("R² = ", round(R2, 3)),
pos = 2)
#Same but by crop Maize
Figure S1
library(dplyr)
library(minpack.lm)
library(nlraa)
# --- Filter for maize only ────────────────────────────────────────────────────
CONV_CA_Ratio_maize <- CONV_CA_Ratio %>% filter(crop == "Maize")
# --- Sort and group ───────────────────────────────────────────────────────────
CONV_CA_Ratio_maize <- CONV_CA_Ratio_maize[order(CONV_CA_Ratio_maize$Y_ratio), ]
CONV_CA_Ratio_maize$group15 <- rep(1:ceiling(nrow(CONV_CA_Ratio_maize)/15), each = 15)[1:nrow(CONV_CA_Ratio_maize)]
grouped_data <- CONV_CA_Ratio_maize %>%
group_by(group15) %>%
summarise(
mean_Y_ratio = mean(Y_ratio, na.rm = TRUE),
mean_DQ = mean(DQ_future_adoption, na.rm = TRUE),
.groups = "drop"
)%>%
filter(mean_Y_ratio <= 2.5) # ← remove groups with mean > 2.5
# --- Fit linear-plateau model ─────────────────────────────────────────────────
fit_lp <- nlsLM(
mean_DQ ~ SSlinp(mean_Y_ratio, a, b, xs),
data = grouped_data,
start = list(a = 0.3, b = 0.5, xs = 1.0),
lower = c(a = 0, b = 0, xs = min(grouped_data$mean_Y_ratio)),
upper = c(a = 1, b = 5, xs = max(grouped_data$mean_Y_ratio))
)
summary(fit_lp)##
## Formula: mean_DQ ~ SSlinp(mean_Y_ratio, a, b, xs)
##
## Parameters:
## Estimate Std. Error t value Pr(>|t|)
## a 0.25105 0.05197 4.831 5.48e-06 ***
## b 0.45801 0.06537 7.007 4.12e-10 ***
## xs 1.19854 0.12223 9.806 6.50e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.1306 on 91 degrees of freedom
##
## Number of iterations to convergence: 7
## Achieved convergence tolerance: 1.49e-08# --- Extract parameters ───────────────────────────────────────────────────────
a_est <- coef(fit_lp)["a"]
b_est <- coef(fit_lp)["b"]
xs_est <- coef(fit_lp)["xs"]
plateau <- a_est + b_est * xs_est
# --- R² ───────────────────────────────────────────────────────────────────────
ss_tot <- sum((grouped_data$mean_DQ - mean(grouped_data$mean_DQ))^2)
R2 <- 1 - sum(residuals(fit_lp)^2) / ss_tot
cat("Breakpoint =", round(xs_est, 3), "\n")## Breakpoint = 1.199cat("Plateau =", round(plateau, 3), "\n")## Plateau = 0.8cat("R² =", round(R2, 3), "\n")## R² = 0.42# --- Target: 15% below plateau ────────────────────────────────────────────────
target <- plateau - 0.15
Y_ratio_target <- (target - a_est) / b_est
cat("Y_ratio for target =", round(Y_ratio_target, 3), "\n")## Y_ratio for target = 0.871# --- Prediction grid ──────────────────────────────────────────────────────────
newdata <- data.frame(
mean_Y_ratio = seq(min(grouped_data$mean_Y_ratio),
max(grouped_data$mean_Y_ratio),
length.out = 500)
)
newdata$pred <- predict(fit_lp, newdata)
# --- Plot ─────────────────────────────────────────────────────────────────────
plot(grouped_data$mean_Y_ratio,
grouped_data$mean_DQ,
pch = 16,
main = "Maize — CA:CT Yield Ratio vs Future Adoption",
xlab = "Mean CA:CT Yield Ratio (groups of 15)",
ylab = "Proportion of future Adoption",
ylim = c(0, 1))
text(x = par("usr")[1], # ← left edge of plot area
y = par("usr")[4], # ← top edge of plot area
labels = "Maize",
adj = c(-0.3, 1.3), # ← small offset from corner
cex = 1.8, font = 2)
lines(newdata$mean_Y_ratio, newdata$pred,
col = "blue", lwd = 2)
# --- Breakpoint lines ─────────────────────────────────────────────────────────
abline(v = xs_est, col = "grey", lwd = 3, lty = 2)
abline(h = plateau, col = "grey", lwd = 2, lty = 3)
points(xs_est, plateau, pch = 19, col = "grey")
# --- Target lines ─────────────────────────────────────────────────────────────
abline(h = target, col = "#009900", lwd = 3, lty = 2)
abline(v = Y_ratio_target, col = "#009900", lwd = 3, lty = 1)
points(Y_ratio_target, target, pch = 19, col = "#009900", cex = 1.5)
# --- Annotations ──────────────────────────────────────────────────────────────
text(xs_est, 0.0,
labels = paste0("Breakpoint = ", round(xs_est, 3)),
col = "black", pos = 4)
text(xs_est, plateau + 0.05,
labels = paste0("Adoption at Plateau = ", round(plateau, 3)),
col = "black", pos = 4)
text(Y_ratio_target, 0.2,
labels = paste0("Yield Ratio (Plateau − 0.15) = ", round(Y_ratio_target, 3)),
col = "#009900", pos = 4)
text(max(grouped_data$mean_Y_ratio), 0.95,
labels = paste0("R² = ", round(R2, 3)),
pos = 2)
Beans
library(dplyr)
library(minpack.lm)
library(nlraa)
# --- Filter for maize only ────────────────────────────────────────────────────
CONV_CA_Ratio_maize <- CONV_CA_Ratio %>% filter(crop == "Beans")
# --- Sort and group ───────────────────────────────────────────────────────────
CONV_CA_Ratio_maize <- CONV_CA_Ratio_maize[order(CONV_CA_Ratio_maize$Y_ratio), ]
CONV_CA_Ratio_maize$group15 <- rep(1:ceiling(nrow(CONV_CA_Ratio_maize)/15), each = 15)[1:nrow(CONV_CA_Ratio_maize)]
grouped_data <- CONV_CA_Ratio_maize %>%
group_by(group15) %>%
summarise(
mean_Y_ratio = mean(Y_ratio, na.rm = TRUE),
mean_DQ = mean(DQ_future_adoption, na.rm = TRUE),
.groups = "drop"
) %>%
filter(mean_Y_ratio <= 2.5) # ← remove groups with mean > 2.5
# --- Fit linear-plateau model ─────────────────────────────────────────────────
fit_lp <- nlsLM(
mean_DQ ~ SSlinp(mean_Y_ratio, a, b, xs),
data = grouped_data,
start = list(a = 0.3, b = 0.5, xs = 1.0),
lower = c(a = 0, b = 0, xs = min(grouped_data$mean_Y_ratio)),
upper = c(a = 1, b = 5, xs = max(grouped_data$mean_Y_ratio))
)
summary(fit_lp)##
## Formula: mean_DQ ~ SSlinp(mean_Y_ratio, a, b, xs)
##
## Parameters:
## Estimate Std. Error t value Pr(>|t|)
## a 0.09695 0.04558 2.127 0.0362 *
## b 0.71400 0.06438 11.090 <2e-16 ***
## xs 1.09049 0.06180 17.646 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.1242 on 90 degrees of freedom
##
## Number of iterations to convergence: 6
## Achieved convergence tolerance: 1.49e-08# --- Extract parameters ───────────────────────────────────────────────────────
a_est <- coef(fit_lp)["a"]
b_est <- coef(fit_lp)["b"]
xs_est <- coef(fit_lp)["xs"]
plateau <- a_est + b_est * xs_est
# --- R² ───────────────────────────────────────────────────────────────────────
ss_tot <- sum((grouped_data$mean_DQ - mean(grouped_data$mean_DQ))^2)
R2 <- 1 - sum(residuals(fit_lp)^2) / ss_tot
cat("Breakpoint =", round(xs_est, 3), "\n")## Breakpoint = 1.09cat("Plateau =", round(plateau, 3), "\n")## Plateau = 0.876cat("R² =", round(R2, 3), "\n")## R² = 0.685# --- Target: 15% below plateau ────────────────────────────────────────────────
target <- plateau - 0.15
Y_ratio_target <- (target - a_est) / b_est
cat("Y_ratio for target =", round(Y_ratio_target, 3), "\n")## Y_ratio for target = 0.88# --- Prediction grid ──────────────────────────────────────────────────────────
newdata <- data.frame(
mean_Y_ratio = seq(min(grouped_data$mean_Y_ratio),
max(grouped_data$mean_Y_ratio),
length.out = 500)
)
newdata$pred <- predict(fit_lp, newdata)
# --- Plot ─────────────────────────────────────────────────────────────────────
plot(grouped_data$mean_Y_ratio,
grouped_data$mean_DQ,
pch = 16,
main = "Beans — CA:CT Yield Ratio vs Future Adoption",
xlab = "Mean CA:CT Yield Ratio (groups of 15)",
ylab = "Proportion of future Adoption",
ylim = c(0, 1))
text(x = par("usr")[1], # ← left edge of plot area
y = par("usr")[4], # ← top edge of plot area
labels = "Beans",
adj = c(-0.3, 1.3), # ← small offset from corner
cex = 1.8, font = 2)
lines(newdata$mean_Y_ratio, newdata$pred,
col = "blue", lwd = 2)
# --- Breakpoint lines ─────────────────────────────────────────────────────────
abline(v = xs_est, col = "grey", lwd = 3, lty = 2)
abline(h = plateau, col = "grey", lwd = 2, lty = 3)
points(xs_est, plateau, pch = 19, col = "grey")
# --- Target lines ─────────────────────────────────────────────────────────────
abline(h = target, col = "#009900", lwd = 3, lty = 2)
abline(v = Y_ratio_target, col = "#009900", lwd = 3, lty = 1)
points(Y_ratio_target, target, pch = 19, col = "#009900", cex = 1.5)
# --- Annotations ──────────────────────────────────────────────────────────────
text(xs_est, 0.0,
labels = paste0("Breakpoint = ", round(xs_est, 3)),
col = "black", pos = 4)
text(xs_est, plateau + 0.05,
labels = paste0("Adoption at Plateau = ", round(plateau, 3)),
col = "black", pos = 4)
text(Y_ratio_target, 0.2,
labels = paste0("Yield Ratio (Plateau − 0.15) = ", round(Y_ratio_target, 3)),
col = "#009900", pos = 4)
text(max(grouped_data$mean_Y_ratio), 0.99,
labels = paste0("R² = ", round(R2, 3)),
pos = 2)
Lets explore also how it looks when we target different slopes levels
## Based on the focus groups discussions
library(dplyr)
library(ggplot2)
# --- Create slope categories ---
CONV_CA_Ratio <- CONV_CA_Ratio %>%
mutate(
slope_cat = case_when(
Slope <= 5 ~ "Low (≤5)",
Slope > 5 & Slope <= 15 ~ "Mid (5–15)",
Slope > 15 ~ "High (>15)"
)
)
# --- Sort within slope category ---
CONV_CA_Ratio <- CONV_CA_Ratio %>%
arrange(slope_cat, Y_ratio)
# --- Create groups of 15 within each slope category ---
CONV_CA_Ratio <- CONV_CA_Ratio %>%
group_by(slope_cat) %>%
mutate(group15 = ceiling(row_number() / 20)) %>%
ungroup()
# --- Aggregate ---
grouped_data <- CONV_CA_Ratio %>%
group_by(slope_cat, group15) %>%
summarise(
mean_Y_ratio = mean(Y_ratio, na.rm = TRUE),
successes = sum(DQ_future_adoption, na.rm = TRUE),
n = n(),
.groups = "drop"
) %>%
mutate(
failures = n - successes,
mean_DQ = successes / n
)
# --- Fit logistic model per slope category ---
models <- grouped_data %>%
group_by(slope_cat) %>%
group_map(~ glm(cbind(successes, failures) ~ mean_Y_ratio,
family = binomial,
data = .x))
names(models) <- unique(grouped_data$slope_cat)
# --- Compute threshold Y_ratio where p = 0.7 for each slope ---
thresholds <- lapply(models, function(m) {
beta0 <- coef(m)[1]
beta1 <- coef(m)[2]
logit_target <- log(0.7 / (1 - 0.7))
(logit_target - beta0) / beta1
})
threshold_df <- data.frame(
slope_cat = names(thresholds),
Y_ratio_07 = unlist(thresholds),
prop = 0.7
)
# --- Plot with facets ---
ggplot(grouped_data, aes(x = mean_Y_ratio, y = mean_DQ)) +
geom_point() +
stat_smooth(method = "glm",
method.args = list(family = "binomial"),
se = FALSE,
color = "blue") +
geom_hline(yintercept = 0.7, color = "red", linetype = "dashed") +
geom_vline(data = threshold_df,
aes(xintercept = Y_ratio_07),
color = "darkgreen",
linewidth = 1) +
geom_point(data = threshold_df,
aes(x = Y_ratio_07, y = prop),
inherit.aes = FALSE,
color = "black") +
geom_text(data = threshold_df,
aes(x = Y_ratio_07, y = 0.15,
label = paste0("Y_ratio = ",
round(Y_ratio_07, 3))),
color = "darkgreen",
inherit.aes = FALSE) +
facet_wrap(~ slope_cat) +
labs(
x = "Mean Y_ratio (groups of 15)",
y = "Proportion DG_future_adoption"
) +
ylim(0,1) +
theme_minimal()
What about the effect of CA on erosion
library(dplyr)
library(tidyr)
# --- Create slope categories and remove NA ---
erosion_table <- CONV_CA_Ratio %>%
filter(
!is.na(DQ_erosion_reduction_rate),
!is.na(Slope),
!is.na(IQR_Season)
) %>%
mutate(
slope_cat = case_when(
Slope <= 5 ~ "Low",
Slope > 5 & Slope <= 15 ~ "Mid",
Slope > 15 ~ "High"
)
) %>%
group_by(IQR_Season, slope_cat) %>%
summarise(
N = n(),
perc = mean(DQ_erosion_reduction_rate) * 100,
.groups = "drop"
) %>%
mutate(
value = paste0(round(perc, 1), "% (", N, ")")
) %>%
dplyr::select(IQR_Season, slope_cat, value) %>%
pivot_wider(
names_from = slope_cat,
values_from = value
) %>%
arrange(IQR_Season)
erosion_tablelibrary(ggplot2)
ggplot(CONV_CA_Ratio, aes(x = IQR_Season, y = Erosion_quanti, group = slope_cat, color = slope_cat)) +
stat_summary(fun = mean, geom = "line", linewidth = 1) +
stat_summary(fun = mean, geom = "point", size = 2) +
labs(
x = "Season",
y = "Erosion_quanti (C - A)",
color = "Slope category",
title = "Evolution of Erosion Reduction Across Seasons by Slope Category"
) +
theme_minimal()
FIgure 3
library(dplyr)
library(ggplot2)
library(forcats)
library(ggh4x) # install.packages("ggh4x") ← only new dependency
plot_df <- CONV_CA_Ratio %>%
filter(!is.na(Erosion_quali), !is.na(slope_cat)) %>%
mutate(
slope_cat = fct_relevel(slope_cat,
"Low (≤5)", "Mid (5–15)", "High (>15)"),
slope_group = "Slope Level (%)", # ← new grouping column
Erosion_quali = fct_relevel(Erosion_quali,
"high_reduction",
"mid_reduction",
"no_reduction"),
Erosion_quali = dplyr::recode(Erosion_quali,
"high_reduction" = "High reduction",
"mid_reduction" = "Moderate reduction",
"no_reduction" = "No reduction")
)
ggplot(plot_df, aes(x = IQR_Season, fill = Erosion_quali)) +
geom_bar(position = "fill") +
facet_nested(~ slope_group + slope_cat) +
scale_y_continuous(labels = scales::percent) +
labs(
x = "Season",
y = "Percentage of Farmers",
fill = "Erosion Reduction\n(Perceived by Farmers)"
) +
theme_minimal(base_size = 14) + # ← raises all text globally
theme(
strip.text = element_text(face = "bold", size = 15), # facet labels
axis.text = element_text(size = 13), # axis tick labels
axis.title = element_text(size = 14), # axis titles
legend.text = element_text(size = 13), # legend items
legend.title = element_text(size = 14), # legend title
panel.spacing = unit(1, "lines")
)
### lets add a new variable of CA performacne classification but with a
threshold of 0.9
library(dplyr)
library(tidyr)
# Step 1: Define target seasons
seasons <- c("23A", "23B", "24A", "24B", "25A", "25B")
# Step 2: Compute number of valid seasons per block
blocks_with_enough_seasons <- CONV_CA_Ratio %>%
filter(IQR_Season %in% seasons, is.finite(Y_ratio)) %>%
group_by(IQR_block) %>%
summarise(n_valid_seasons = n_distinct(IQR_Season), .groups = "drop") %>%
filter(n_valid_seasons >= 4)
# Step 3: Filter main dataset to include only blocks with ≥4 seasons
filtered_yield <- CONV_CA_Ratio %>%
filter(IQR_block %in% blocks_with_enough_seasons$IQR_block,
IQR_Season %in% seasons,
is.finite(Y_ratio))
# Step 4: Prepare wide-format Relative Yield table
RY_wide <- filtered_yield %>%
group_by(IQR_block, IQR_Season) %>%
summarise(Y_ratio = mean(Y_ratio, na.rm = TRUE), .groups = "drop") %>%
pivot_wider(
names_from = IQR_Season,
values_from = Y_ratio,
names_prefix = "Y_ratio_"
)
# Step 5: Compute performance typology
Perf_Type_09 <- filtered_yield %>%
group_by(IQR_block) %>%
summarise(
n_seasons = n_distinct(IQR_Season),
any_below_0.9 = any(Y_ratio < 0.9),
mean_ratio = mean(Y_ratio, na.rm = TRUE),
.groups = "drop"
) %>%
mutate(
CA_typology_09 = case_when(
!any_below_0.9 ~ "Good",
any_below_0.9 & mean_ratio >= 0.9 ~ "Acceptable",
any_below_0.9 & mean_ratio < 0.9 & mean_ratio > 0.8 ~ "Poor",
mean_ratio <= 0.8 ~ "Very Poor"
)
)#Figure 4B
library(ggplot2)
library(forcats)
# Collapse Good into Acceptable and set order
Perf_Type_09 <- Perf_Type_09 %>%
mutate(CA_typology_09 = fct_collapse(CA_typology_09,
Acceptable = c("Acceptable", "Good")),
CA_typology_09 = factor(CA_typology_09,
levels = c("Very Poor", "Poor", "Acceptable")))
# Overall percent bar
ggplot(Perf_Type_09, aes(x = CA_typology_09, fill = CA_typology_09)) +
geom_bar(aes(y = after_stat(count / sum(count) * 100))) +
# ── Panel label B ──────────────────────────────────────────────────────────
annotate("text",
x = Inf, y = Inf, # ← right edge
hjust = 1.3, vjust = 1.5, # ← offset from top right corner
label = "B",
size = 7, fontface = "bold",
colour = "black") +
labs(
x = "Performance along six crop seasons",
y = "Percentage of Sites"
) +
scale_fill_manual(values = c(
"Very Poor" = "#CC0000",
"Poor" = "#FF9900",
"Acceptable" = "#009900"
)) +
theme_minimal(base_size = 16) +
theme(
legend.position = "none",
panel.grid.major = element_blank(), # ← removes major grid lines
panel.grid.minor = element_blank(), # ← removes minor grid lines
axis.line = element_line(colour = "black"), # ← adds axis lines
axis.ticks = element_line(colour = "black") # ← adds tick marks
)
Figure 2 with overall
library(dplyr)
library(ggplot2)
library(glue)
# 1. Mean Y_ratio per block x crop ────────────────────────────────────────────
block_means_crop <- CONV_CA_Ratio %>%
filter(is.finite(Y_ratio), !is.na(IQF_agzone), !is.na(crop)) %>%
group_by(IQR_block, IQF_agzone, crop) %>%
summarise(mean_Y_ratio = mean(Y_ratio, na.rm = TRUE), .groups = "drop") %>%
mutate(IQF_agzone = droplevels(as.factor(IQF_agzone)))
# 2. Overall (across crops) ───────────────────────────────────────────────────
block_means_overall <- CONV_CA_Ratio %>%
filter(is.finite(Y_ratio), !is.na(IQF_agzone)) %>%
group_by(IQR_block, IQF_agzone) %>%
summarise(mean_Y_ratio = mean(Y_ratio, na.rm = TRUE), .groups = "drop") %>%
mutate(crop = "Overall",
IQF_agzone = droplevels(as.factor(IQF_agzone)))
# 3. Combine ───────────────────────────────────────────────────────────────────
block_means <- bind_rows(block_means_crop, block_means_overall) %>%
mutate(crop = factor(crop, levels = c("Maize", "Beans", "Overall")))
# 4. Order agzones by overall mean yield ratio ────────────────────────────────
agzone_order <- block_means_overall %>%
group_by(IQF_agzone) %>%
summarise(agzone_mean = mean(mean_Y_ratio, na.rm = TRUE), .groups = "drop") %>%
arrange(agzone_mean) %>%
pull(IQF_agzone)
block_means <- block_means %>%
mutate(IQF_agzone = factor(IQF_agzone, levels = agzone_order))
# 4b. Append "All zones" (pools all zones, displayed last) ────────────────────
all_zones <- block_means %>%
mutate(IQF_agzone = "All zones")
block_means <- bind_rows(block_means, all_zones) %>%
mutate(IQF_agzone = factor(IQF_agzone,
levels = c(as.character(agzone_order), "All zones")))
# 5. Summary statistics per agzone x crop ─────────────────────────────────────
agzone_summary <- block_means %>%
group_by(IQF_agzone, crop) %>%
summarise(
avg_ratio = mean(mean_Y_ratio, na.rm = TRUE),
n_total = n(),
n_ge1 = sum(mean_Y_ratio >= 1, na.rm = TRUE),
percent_ge1 = round(100 * n_ge1 / n_total),
.groups = "drop"
) %>%
mutate(
label_top = paste0(" ", percent_ge1),
label_bottom = paste0(" ", formatC(avg_ratio, format = "f", digits = 2)),
y_top = 0.15,
y_bottom = 0.01
)
# 6. Header labels (≥1% and μ) per facet ──────────────────────────────────────
header_labels <- data.frame(
crop = factor(c("Maize", "Beans", "Overall"),
levels = c("Maize", "Beans", "Overall")),
x = 0.7,
y_top = 0.15,
y_bottom = 0.01
)
# 7. Bold face vector for x-axis ("All zones" in bold, others plain) ──────────
agzone_levels <- levels(block_means$IQF_agzone)
axis_face <- ifelse(agzone_levels == "All zones", "bold", "plain")
# Create a named vector for the labeller
crop_labels <- c("Maize" = "A) Maize", "Beans" = "B) Beans", "Overall" = "C) Overall")
# 8. Plot ──────────────────────────────────────────────────────────────────────
ggplot(block_means, aes(x = IQF_agzone, y = mean_Y_ratio)) +
geom_boxplot(
outlier.shape = NA,
fill = "gray90",
color = "black",
width = 0.5
) +
geom_hline(
yintercept = 1,
linetype = "dashed",
color = "black",
linewidth = 0.6
) +
geom_jitter(
width = 0.2,
alpha = 0.4,
size = 1.2,
color = "grey60"
) +
geom_text(
data = agzone_summary,
aes(x = IQF_agzone, y = y_top, label = label_top),
size = 5
) +
geom_text(
data = agzone_summary,
aes(x = IQF_agzone, y = y_bottom, label = label_bottom),
size = 5
) +
geom_text(
data = header_labels,
aes(x = x, y = y_top, label = "≥1(%)"),
hjust = 1, size = 5, inherit.aes = FALSE
) +
geom_text(
data = header_labels,
aes(x = x, y = y_bottom, label = "μ"),
hjust = 1, size = 5, inherit.aes = FALSE
) +
# Vertical separator line before "All zones"
geom_vline(
xintercept = length(agzone_levels) - 0.5,
linetype = "dotted",
color = "black",
linewidth = 0.5
) +
facet_wrap(~ crop, ncol = 1, labeller = labeller(crop = crop_labels)) +
coord_cartesian(ylim = c(0, 1.65), clip = "off") +
scale_y_continuous(breaks = seq(0, 1.6, by = 0.2)) +
scale_x_discrete(expand = expansion(mult = c(0.12, 0.02))) +
labs(
x = "Agroecological Zone",
y = "Mean Yield Ratio (CA:CT)"
) +
theme_minimal(base_size = 20) +
theme(
strip.text = element_text(face = "bold", size = 20),
strip.background = element_rect(fill = "grey90", colour = NA),
axis.title.x = element_text(size = 20),
axis.title.y = element_text(size = 20),
axis.text.x = element_text(angle = 45, hjust = 1, size = 20,
face = axis_face),
axis.text.y = element_text(size = 20),
axis.line = element_line(color = "black"),
axis.ticks = element_line(color = "black"),
panel.grid.major = element_blank(),
panel.grid.minor = element_blank(),
plot.margin = margin(20, 20, 20, 20)
)
library(dplyr)
library(ggplot2)
library(glue)
# 1. Mean Y_ratio per block (across seasons) and remove NA districts
block_means <- CONV_CA_Ratio %>%
filter(is.finite(Y_ratio), !is.na(IQF_District)) %>%
group_by(IQR_block, IQF_District, IQF_environment) %>%
summarise(mean_Y_ratio = mean(Y_ratio, na.rm = TRUE), .groups = "drop") %>%
mutate(IQF_District = droplevels(as.factor(IQF_District)))
# 2. Order districts by mean yield ratio
district_order <- block_means %>%
group_by(IQF_District) %>%
summarise(district_mean = mean(mean_Y_ratio, na.rm = TRUE), .groups = "drop") %>%
arrange(district_mean) %>%
pull(IQF_District)
block_means <- block_means %>%
mutate(IQF_District = factor(IQF_District, levels = district_order))
# 3. Summary statistics
district_summary <- block_means %>%
group_by(IQF_District) %>%
summarise(
avg_ratio = mean(mean_Y_ratio, na.rm = TRUE),
n_total = n(),
n_ge1 = sum(mean_Y_ratio >= 1, na.rm = TRUE),
percent_ge1 = round(100 * n_ge1 / n_total),
.groups = "drop"
) %>%
mutate(
label_top = paste0(" ", percent_ge1),
label_bottom = paste0(" ", formatC(avg_ratio, format = "f", digits = 2)),
y_top = 0.08,
y_bottom = 0.02
)
# 4. Plot
ggplot(block_means, aes(x = IQF_District, y = mean_Y_ratio)) +
geom_boxplot(
outlier.shape = NA,
fill = "gray90",
color = "black",
width = 0.5
) +
geom_jitter(
aes(color = IQF_environment),
width = 0.2,
alpha = 0.6,
size = 2
) +
geom_text(
data = district_summary,
aes(x = IQF_District, y = y_top, label = label_top),
size = 3
) +
geom_text(
data = district_summary,
aes(x = IQF_District, y = y_bottom, label = label_bottom),
size = 3
) +
coord_cartesian(ylim = c(0, 1.65), clip = "off") +
scale_y_continuous(
breaks = seq(0, 1.6, by = 0.2)
) +
scale_x_discrete(
expand = expansion(mult = c(0.12, 0.02))
) +
labs(
x = "District",
y = "Mean Yield Ratio (CA:CT)",
color = "Environment"
) +
theme_minimal(base_size = 12) +
theme(
axis.title.x = element_text(size = 14),
axis.title.y = element_text(size = 14),
axis.text.x = element_text(angle = 45, hjust = 1),
axis.line = element_line(color = "black"),
axis.ticks = element_line(color = "black"),
panel.grid.major = element_blank(),
panel.grid.minor = element_blank(),
plot.margin = margin(20, 20, 20, 20),
legend.title = element_text(size = 12),
legend.text = element_text(size = 10),
legend.position = "right",
legend.box.margin = margin(0,0,0,0),
legend.spacing.y = unit(0.2, "cm"),
legend.key.height = unit(0.35, "cm")
) +
annotate("text", x = 0.3, y = 0.08, label = "≥1(%)", hjust = 1, size = 3) +
annotate("text", x = 0.3, y = 0.02, label = "μ", hjust = 1, size = 3)
CA performance on individual plots rather than averages
Description of the analysis: We analyzed the block-level average performance of Conservation Agriculture (CA) relative to Conventional Agriculture (CONV) by calculating the mean CA:CONV yield ratio per block across all seasons (i.e., 3 years = 6 seasons, 3 maize + 3 beans). A cumulative distribution plot was used to visualize how many blocks reached or exceeded yield parity (Y_ratio ≥ 1). Additionally, density plots by agro-environment visualized how CA performance varied across zones.
Description of the results: As we showed previously, when examining the average performance of CA at the country, environment, or district level, results are discouraging, suggesting that a blanket recommendation of this package is unlikely to generate broad impact or adoption, even if we target the best environmetn or district. However, when zooming in to the plot level, there are a substantial number of cases where CA either outperforms CONV or has a similar performance. Specifically, CA outperforms CONV in 18% of blocks, and in 21% of cases achieves a yield ratio close to or above acceptable thresholds (Shown before in the figure of Performance Typologies)
library(dplyr)
library(ggplot2)
# 1. Summarize data to block level
block_data <- CONV_CA_Ratio %>%
filter(is.finite(Y_ratio)) %>%
group_by(IQR_block) %>%
summarise(
avg_Y_ratio = mean(Y_ratio, na.rm = TRUE),
IQF_environment = first(IQF_environment), # to retain for optional grouping
.groups = "drop"
) %>%
arrange(avg_Y_ratio) %>%
mutate(cum_percent = 100 * (row_number() / n()))
# 2. Plot cumulative distribution
ggplot(block_data, aes(x = avg_Y_ratio, y = cum_percent)) +
geom_point(size = 2, alpha = 0.7, color = "darkorange") +
geom_vline(xintercept = 1, linetype = "dashed", color = "blue") +
scale_x_continuous(limits = c(0, 2.5)) +
scale_y_continuous(breaks = seq(0, 100, by = 10)) +
labs(
x = "Block-Level Average Yield Ratio (CA / CONV)",
y = "Cumulative % of Blocks",
title = "Cumulative Distribution of Relative Yield (CA vs CONV)",
subtitle = "Each point is a block's average over all seasons; Dashed line = yield parity",
) +
theme_minimal() +
theme(
plot.title = element_text(face = "bold", hjust = 0.5),
plot.subtitle = element_text(size = 10, hjust = 0.5)
)
##Visual distribution of CONV advantage vs CA_Yield
library(ggplot2)
ggplot(CONV_CA_Ratio, aes(x = Y_ratio, fill = IQF_environment)) +
geom_density(alpha = 0.3) +
geom_vline(xintercept = 1, linetype = "dashed") +
labs(title = "Density of Relative Yield (CA / CONV) by Agzone")
Drivers of yield performance
#Explore, visually the protential drivers and the relatio nwith Y_ration of the blocks
Cathegoric variables
Variables we have so far (after filtering some out and re-grouping cathegories with low N: ICR_Org_C_avg num: ICR_K_perc_23A num: ICR_Avail_P_avg num: ICR_arable_land_avg num: ICF_Alt_avg num: ICR_rainfall_avg num: Comp_amount_corr_quali num: ICF_planting_date num: IQR_TF_tubura_client Factor: IQR_plot_fert_quality Factor: IQR_soil_texture Factor w/ 4 levels: IQR_soil_color Factor w/ 6 levels: IQF_position_slope Factor w/ 7 levels: Weed_species_A_combined: Factor w/ 5 levels: IQF_environment chr “Mid_Low_Wet” “Mid_Low_Wet” “Mid_Low_Wet” “Mid_Low_Wet” Output variables: avg_Y_ratio num: CA_typology factor:
CA Performance as affected by ICR_Org_C_avg
with avg_Y_ratio
Notes: when we look at all the database pooled together, there is no relation at all among SOC and performance of the CA. Howeever, where we split the analysis be Environment, then an interation becomes clear, where there is no relation in mid altitude and the relation is quite significant and positive in Highland and sligtly negative in transition. Conclusion. SOC should remain for the drivers analysis
library(ggplot2)
# 1. Fit linear regression model
model <- lm(avg_Y_ratio ~ ICR_Org_C_avg, data = drivers_data)
# 2. Extract p-value from model summary
p_val <- summary(model)$coefficients[2, 4] # p-value for slope
# 3. Create plot
ggplot(drivers_data, aes(x = ICR_Org_C_avg, y = avg_Y_ratio)) +
geom_point(alpha = 0.6, color = "steelblue") +
geom_smooth(method = "lm", se = TRUE, color = "darkred", linetype = "dashed") +
labs(
title = "ICR_Org_C_avg vs. avg_Y_ratio",
x = "Organic Carbon (%)",
y = "Average Yield Ratio"
) +
annotate("text",
x = Inf, y = Inf, hjust = 1.1, vjust = 1.5,
label = paste0("p = ", signif(p_val, 3)),
size = 5, color = "black") +
theme_minimal()
### What if we facet by environment? to see if interactions are masking
the effect
library(ggplot2)
library(dplyr)
library(purrr)
library(tidyr)
#----------------------------------------------
# 1. Fit separate linear models per environment
#----------------------------------------------
stats_df <- drivers_data %>%
filter(!is.na(ICR_Org_C_avg), !is.na(avg_Y_ratio), !is.na(IQF_environment)) %>%
group_by(IQF_environment) %>%
nest() %>%
mutate(
model = map(data, ~ lm(avg_Y_ratio ~ ICR_Org_C_avg, data = .x)),
p_value = map_dbl(model, ~ summary(.x)$coefficients[2, 4]),
p_label = paste0("p = ", signif(p_value, 3))
) %>%
dplyr::select(IQF_environment, p_label)
#----------------------------------------------
# 2. Prepare annotation positions (adjust if needed)
#----------------------------------------------
# Choose fixed position for annotation
label_y <- max(drivers_data$avg_Y_ratio, na.rm = TRUE) * 0.95
label_x <- max(drivers_data$ICR_Org_C_avg, na.rm = TRUE) * 0.95
stats_df <- stats_df %>%
mutate(x = label_x, y = label_y)
#----------------------------------------------
# 3. Plot: points, lm line, facet, p-values
#----------------------------------------------
ggplot(drivers_data, aes(x = ICR_Org_C_avg, y = avg_Y_ratio)) +
geom_point(alpha = 0.6, color = "steelblue") +
geom_smooth(method = "lm", se = TRUE, color = "darkred", linetype = "dashed") +
geom_text(
data = stats_df,
aes(x = x, y = y, label = p_label),
inherit.aes = FALSE,
size = 5,
hjust = 1,
vjust = 1
) +
facet_wrap(~ IQF_environment) +
labs(
title = "ICR_Org_C_avg vs. avg_Y_ratio",
subtitle = "Linear regressions and p-values shown per environment",
x = "Organic Carbon (%)",
y = "Average Yield Ratio"
) +
theme_minimal()
with CA_typology
library(ggplot2)
library(dplyr)
library(multcompView)
# 1. Fit ANOVA model
model_anova <- aov(ICR_Org_C_avg ~ CA_typology, data = drivers_data)
# 2. Tukey HSD post-hoc test
tukey <- TukeyHSD(model_anova)
# 3. Extract group letters
letters_df <- multcompLetters(tukey$CA_typology[,"p adj"])$Letters
letters_df <- data.frame(CA_typology = names(letters_df), Letters = letters_df)
# 4. Calculate group means
means_df <- drivers_data %>%
group_by(CA_typology) %>%
summarise(mean_val = mean(ICR_Org_C_avg, na.rm = TRUE), .groups = "drop") %>%
mutate(mean_label = paste0("mean = ", round(mean_val, 2)))
# 5. Combine letters and means
plot_labels <- left_join(letters_df, means_df, by = "CA_typology")
# 6. Create boxplot with letters and means at fixed Y positions
ggplot(drivers_data, aes(x = CA_typology, y = ICR_Org_C_avg)) +
geom_boxplot(outlier.shape = NA, fill = "lightgray", color = "black") +
geom_jitter(width = 0.2, alpha = 0.6, color = "steelblue", size = 1.5) +
geom_text(data = plot_labels, aes(x = CA_typology, y = 7.5, label = Letters),
fontface = "bold", size = 5, color = "darkred") +
geom_text(data = plot_labels, aes(x = CA_typology, y = 7, label = mean_label),
fontface = "plain", size = 4, color = "black") +
labs(
title = "ICR_Org_C_avg by CA_typology",
subtitle = "Letters indicate significant differences (Tukey HSD)",
x = "CA Typology",
y = "Organic Carbon (%)"
) +
theme_minimal()
### What if we facet by Environment? myve some interactions masking the
effect?
library(ggplot2)
library(dplyr)
library(multcompView)
library(tidyr)
library(purrr)
#---------------------------------------------------------
# 1. Function to run ANOVA + Tukey + extract letters
#---------------------------------------------------------
get_letters <- function(df) {
mod <- aov(ICR_Org_C_avg ~ CA_typology, data = df)
tuk <- TukeyHSD(mod)
lets <- multcompLetters(tuk$CA_typology[, "p adj"])$Letters
data.frame(
CA_typology = names(lets),
Letters = lets,
stringsAsFactors = FALSE
)
}
#---------------------------------------------------------
# 2. Run statistics PER ENVIRONMENT
#---------------------------------------------------------
stats_df <- drivers_data %>%
filter(!is.na(CA_typology)) %>%
group_by(IQF_environment) %>%
nest() %>%
mutate(
letters = map(data, get_letters),
means = map(data, ~ .x %>%
group_by(CA_typology) %>%
summarise(mean_val = mean(ICR_Org_C_avg, na.rm = TRUE),
.groups = "drop"))
) %>%
unnest(letters) %>%
left_join(
unnest(dplyr::select(., IQF_environment, means)),
by = c("IQF_environment", "CA_typology")
) %>%
mutate(mean_label = paste0("mean = ", round(mean_val, 2)))
#---------------------------------------------------------
# 3. Plot (letters + means fixed at Y)
#---------------------------------------------------------
ggplot(drivers_data, aes(x = CA_typology, y = ICR_Org_C_avg)) +
geom_boxplot(outlier.shape = NA, fill = "lightgray") +
geom_jitter(width = 0.2, alpha = 0.6, color = "steelblue", size = 1.4) +
# Letters
geom_text(
data = stats_df,
aes(x = CA_typology, y = 9.5, label = Letters),
inherit.aes = FALSE,
fontface = "bold",
size = 3,
color = "darkred"
) +
# Means
geom_text(
data = stats_df,
aes(x = CA_typology, y = 8.25, label = mean_label),
inherit.aes = FALSE,
size = 3
) +
facet_wrap(~ IQF_environment) +
labs(
title = "ICR_Org_C_avg by CA_typology (per environment)",
subtitle = "Letters = Tukey HSD within each environment",
x = "CA Typology",
y = "Organic Carbon (%)"
) +
theme_minimal()
## CA Performance as affected by ICR_K_perc_23A ### with avg_Y_ratio
Notes: there is some minimal relation in highland. When we look at the
correlation among SOC and K, its relatively high in highland. So we
could just use SOC in the drivers analysis and exclude K Conclusion:
Keep SCO but explcude K
library(ggplot2)
library(dplyr)
library(purrr)
library(tidyr)
#----------------------------------------------
# 1. Fit separate linear models per environment
#----------------------------------------------
stats_df <- drivers_data %>%
filter(!is.na(ICR_K_perc_23A), !is.na(avg_Y_ratio), !is.na(IQF_environment)) %>%
group_by(IQF_environment) %>%
nest() %>%
mutate(
model = map(data, ~ lm(avg_Y_ratio ~ ICR_K_perc_23A, data = .x)),
p_value = map_dbl(model, ~ summary(.x)$coefficients[2, 4]),
p_label = paste0("p = ", signif(p_value, 3))
) %>%
dplyr::select(IQF_environment, p_label)
#----------------------------------------------
# 2. Prepare annotation positions
#----------------------------------------------
label_y <- max(drivers_data$avg_Y_ratio, na.rm = TRUE) * 0.95
label_x <- max(drivers_data$ICR_K_perc_23A, na.rm = TRUE) * 0.95
stats_df <- stats_df %>%
mutate(x = label_x, y = label_y)
#----------------------------------------------
# 3. Plot: points, regression line, facet, p-values
#----------------------------------------------
ggplot(drivers_data, aes(x = ICR_K_perc_23A, y = avg_Y_ratio)) +
geom_point(alpha = 0.6, color = "steelblue") +
geom_smooth(method = "lm", se = TRUE, color = "darkred", linetype = "dashed") +
geom_text(
data = stats_df,
aes(x = x, y = y, label = p_label),
inherit.aes = FALSE,
size = 5,
hjust = 1,
vjust = 1
) +
facet_wrap(~ IQF_environment) +
labs(
title = "ICR_K_perc_23A vs. avg_Y_ratio",
subtitle = "Linear regressions and p-values shown per environment",
x = "Potassium (%) in 23A",
y = "Average Yield Ratio"
) +
theme_minimal()
### with CA_typology
library(ggplot2)
library(dplyr)
library(multcompView)
library(tidyr)
library(purrr)
#---------------------------------------------------------
# 1. Function to run ANOVA + Tukey + extract letters
#---------------------------------------------------------
get_letters <- function(df) {
mod <- aov(ICR_K_perc_23A ~ CA_typology, data = df)
tuk <- TukeyHSD(mod)
lets <- multcompLetters(tuk$CA_typology[, "p adj"])$Letters
data.frame(
CA_typology = names(lets),
Letters = lets,
stringsAsFactors = FALSE
)
}
#---------------------------------------------------------
# 2. Run statistics PER ENVIRONMENT
#---------------------------------------------------------
stats_df <- drivers_data %>%
filter(!is.na(CA_typology), !is.na(ICR_K_perc_23A)) %>%
group_by(IQF_environment) %>%
nest() %>%
mutate(
letters = map(data, get_letters),
means = map(data, ~ .x %>%
group_by(CA_typology) %>%
summarise(mean_val = mean(ICR_K_perc_23A, na.rm = TRUE),
.groups = "drop"))
) %>%
unnest(letters) %>%
left_join(
unnest(dplyr::select(., IQF_environment, means)),
by = c("IQF_environment", "CA_typology")
) %>%
mutate(mean_label = paste0("mean = ", round(mean_val, 2)))
#---------------------------------------------------------
# 3. Plot (letters + means fixed at Y)
#---------------------------------------------------------
ggplot(drivers_data, aes(x = CA_typology, y = ICR_K_perc_23A)) +
geom_boxplot(outlier.shape = NA, fill = "lightgray") +
geom_jitter(width = 0.2, alpha = 0.6, color = "steelblue", size = 1.4) +
# Letters
geom_text(
data = stats_df,
aes(x = CA_typology, y = 3.2, label = Letters),
inherit.aes = FALSE,
fontface = "bold",
size = 3,
color = "darkred"
) +
# Means
geom_text(
data = stats_df,
aes(x = CA_typology, y = 2.8, label = mean_label),
inherit.aes = FALSE,
size = 3
) +
facet_wrap(~ IQF_environment) +
labs(
title = "ICR_K_perc_23A by CA_typology (per environment)",
subtitle = "Letters = Tukey HSD within each environment",
x = "CA Typology",
y = "Potassium (%) in 23A"
) +
theme_minimal()
## Since SOC semms to have some explanatory power and K very little,
lets see the correlation among them Notes: there is quite a correlation
among K and SOC in highland
library(ggplot2)
library(dplyr)
library(purrr)
library(tidyr)
#---------------------------------------------------------
# 1. Fit linear models per environment
#---------------------------------------------------------
correlation_stats <- drivers_data %>%
filter(!is.na(ICR_Org_C_avg), !is.na(ICR_K_perc_23A), !is.na(IQF_environment)) %>%
group_by(IQF_environment) %>%
nest() %>%
mutate(
model = map(data, ~ lm(ICR_K_perc_23A ~ ICR_Org_C_avg, data = .x)),
p_value = map_dbl(model, ~ summary(.x)$coefficients[2, 4]),
r_squared = map_dbl(model, ~ summary(.x)$r.squared),
stat_label = paste0("p = ", signif(p_value, 3),
"\nR", "\u00B2", " = ", round(r_squared, 2))
)
#---------------------------------------------------------
# 2. Set annotation position
#---------------------------------------------------------
label_y <- max(drivers_data$ICR_K_perc_23A, na.rm = TRUE) * 0.95
label_x <- max(drivers_data$ICR_Org_C_avg, na.rm = TRUE) * 0.95
correlation_stats <- correlation_stats %>%
mutate(x = label_x, y = label_y)
#---------------------------------------------------------
# 3. Plot: scatter + regression + p + R² + facet
#---------------------------------------------------------
ggplot(drivers_data, aes(x = ICR_Org_C_avg, y = ICR_K_perc_23A)) +
geom_point(alpha = 0.6, color = "steelblue", size = 1.5) +
geom_smooth(method = "lm", se = TRUE, color = "darkred", linetype = "dashed") +
geom_text(
data = correlation_stats,
aes(x = x, y = y, label = stat_label),
inherit.aes = FALSE,
hjust = 1, vjust = 1,
size = 4
) +
facet_wrap(~ IQF_environment) +
labs(
title = "Correlation between Organic Carbon and Potassium (23A)",
subtitle = "Linear regression per environment with p-values and R\u00B2",
x = "ICR_Org_C_avg (Organic Carbon %)",
y = "ICR_K_perc_23A (Potassium %)"
) +
theme_minimal()
### What if instead of contineous, we define deffitient or non
defficiennt
library(dplyr)
library(ggplot2)
library(segmented)
# ===============================
# 1) Average env_index per block across seasons
# ===============================
env_block_avg <- CONV_CA_Ratio %>%
filter(!is.na(IQR_block), !is.na(env_index)) %>%
group_by(IQR_block) %>%
summarise(
env_index_avg = mean(env_index, na.rm = TRUE),
.groups = "drop"
)
# ===============================
# 2) Ensure drivers_data is block-level for K
# ===============================
k_block <- drivers_data %>%
filter(!is.na(IQR_block), !is.na(ICR_K_perc_23A)) %>%
group_by(IQR_block) %>%
summarise(
ICR_K_perc_23A = mean(ICR_K_perc_23A, na.rm = TRUE),
.groups = "drop"
)
# ===============================
# 3) Join datasets
# ===============================
plot_data <- env_block_avg %>%
left_join(k_block, by = "IQR_block") %>%
filter(!is.na(ICR_K_perc_23A))
cat("Blocks used:", nrow(plot_data), "\n")## Blocks used: 459# ===============================
# 4) Fit linear model
# ===============================
lm_model <- lm(env_index_avg ~ ICR_K_perc_23A,
data = plot_data)
# ===============================
# 5) Fit segmented (breakpoint) model
# ===============================
seg_model <- segmented(
lm_model,
seg.Z = ~ ICR_K_perc_23A,
psi = list(ICR_K_perc_23A = median(plot_data$ICR_K_perc_23A, na.rm = TRUE))
)
summary(seg_model)##
## ***Regression Model with Segmented Relationship(s)***
##
## Call:
## segmented.lm(obj = lm_model, seg.Z = ~ICR_K_perc_23A, psi = list(ICR_K_perc_23A = median(plot_data$ICR_K_perc_23A,
## na.rm = TRUE)))
##
## Estimated Break-Point(s):
## Est. St.Err
## psi1.ICR_K_perc_23A 0.3 0.045
##
## Coefficients of the linear terms:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 0.71785 0.06459 11.115 < 2e-16 ***
## ICR_K_perc_23A 1.06639 0.31448 3.391 0.000757 ***
## U1.ICR_K_perc_23A -1.00761 0.31623 -3.186 NA
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.2824 on 455 degrees of freedom
## Multiple R-Squared: 0.08489, Adjusted R-squared: 0.07885
##
## Boot restarting based on 6 samples. Last fit:
## Convergence attained in 1 iterations (rel. change 1.5678e-07)# Extract breakpoint
breakpoint <- seg_model$psi[2]
cat("Estimated Breakpoint:", round(breakpoint, 3), "\n")## Estimated Breakpoint: 0.3# Add predictions
plot_data$predicted <- predict(seg_model)
# ===============================
# 6) Plot
# ===============================
ggplot(plot_data,
aes(x = ICR_K_perc_23A,
y = env_index_avg)) +
geom_point(alpha = 0.7,
color = "steelblue",
size = 2) +
geom_line(aes(y = predicted),
color = "red",
size = 1.2) +
geom_vline(xintercept = breakpoint,
linetype = "dashed",
color = "darkgreen",
size = 1) +
labs(
title = "Breakpoint Model: Environmental Index vs Potassium (%)",
subtitle = paste0("Estimated breakpoint ≈ ",
round(breakpoint, 3)),
x = "ICR_K_perc_23A",
y = "Mean env_index across seasons"
) +
theme_minimal(base_size = 14)
Based on litterature: VLow: <0.2 meq/100g Low: >=0.2 and < 0.5
meq/100g. Opt: =>0.5meq/100g. https://link.springer.com/article/10.1186/s40066-018-0165-5#:~:text=Exchangeable%20K,-The%20soil%20exchangeable&text=The%20range%20of%20values%20recorded,their%20soil%20exchangeable%20K%2C%20respectively.
Then adjusted to our findings: VLow: <0.3 meq/100g Low: >=0.23 and < 0.56 meq/100g. Opt: =>0.6meq/100g.
FIgure of K_Factor with Y_ratio
library(ggplot2)
library(dplyr)
# 1️⃣ Filter valid data
df_filt <- drivers_data %>%
filter(
!is.na(avg_Y_ratio),
!is.na(K_factor)
) %>%
droplevels()
# 2️⃣ Run ANOVA
anova_model <- aov(avg_Y_ratio ~ K_factor, data = df_filt)
anova_p <- summary(anova_model)[[1]][["Pr(>F)"]][1]
# 3️⃣ Run Kruskal-Wallis test
kruskal_p <- kruskal.test(avg_Y_ratio ~ K_factor, data = df_filt)$p.value
# 4️⃣ Build annotation label
test_label <- paste0(
"ANOVA p = ", signif(anova_p, 3),
"\nKruskal p = ", signif(kruskal_p, 3)
)
# 5️⃣ Calculate group means
mean_df <- df_filt %>%
group_by(K_factor) %>%
summarise(
mean_val = mean(avg_Y_ratio, na.rm = TRUE),
.groups = "drop"
) %>%
mutate(mean_label = paste0("mean = ", round(mean_val, 2)))
# 6️⃣ Plot
ggplot(df_filt, aes(x = K_factor, y = avg_Y_ratio)) +
geom_boxplot(outlier.shape = NA, fill = "lightgray") +
geom_jitter(width = 0.2, alpha = 0.6, color = "steelblue") +
# Means under each box
geom_text(
data = mean_df,
aes(
x = K_factor,
y = min(df_filt$avg_Y_ratio, na.rm = TRUE) * 0.95,
label = mean_label
),
inherit.aes = FALSE,
size = 4
) +
# P-values annotation
annotate(
"text",
x = Inf, y = Inf,
hjust = 1.1, vjust = 1.5,
label = test_label,
size = 5
) +
labs(
title = "avg_Y_ratio by K_factor",
subtitle = "ANOVA and Kruskal–Wallis test results shown",
x = "Potassium Level (K_factor)",
y = "Average Yield Ratio"
) +
theme_minimal()
CA Performance as affected by ICR_Avail_P_avg
with avg_Y_ratio
NOtes: the values sees high in general and not relation with the yield of the plots (environmental index calculated as the ) Conclusions: poor
library(ggplot2)
library(dplyr)
library(purrr)
library(tidyr)
#----------------------------------------------
# 1. Fit separate linear models per environment
#----------------------------------------------
stats_df <- drivers_data %>%
filter(!is.na(ICR_Avail_P_avg), !is.na(avg_Y_ratio), !is.na(IQF_environment)) %>%
group_by(IQF_environment) %>%
nest() %>%
mutate(
model = map(data, ~ lm(avg_Y_ratio ~ ICR_Avail_P_avg, data = .x)),
p_value = map_dbl(model, ~ summary(.x)$coefficients[2, 4]),
p_label = paste0("p = ", signif(p_value, 3))
) %>%
dplyr::select(IQF_environment, p_label)
#----------------------------------------------
# 2. Prepare annotation positions
#----------------------------------------------
label_y <- max(drivers_data$avg_Y_ratio, na.rm = TRUE) * 0.95
label_x <- max(drivers_data$ICR_Avail_P_avg, na.rm = TRUE) * 0.95
stats_df <- stats_df %>%
mutate(x = label_x, y = label_y)
#----------------------------------------------
# 3. Plot: points, regression line, facet, p-values
#----------------------------------------------
ggplot(drivers_data, aes(x = ICR_Avail_P_avg, y = avg_Y_ratio)) +
geom_point(alpha = 0.6, color = "steelblue") +
geom_smooth(method = "lm", se = TRUE, color = "darkred", linetype = "dashed") +
geom_text(
data = stats_df,
aes(x = x, y = y, label = p_label),
inherit.aes = FALSE,
size = 5,
hjust = 1,
vjust = 1
) +
facet_wrap(~ IQF_environment) +
labs(
title = "ICR_Avail_P_avg vs. avg_Y_ratio",
subtitle = "Linear regressions and p-values shown per environment",
x = "Available Phosphorus (ppm)",
y = "Average Yield Ratio"
) +
theme_minimal()
###ICR_Avail_P_avg by CA_typology
library(ggplot2)
library(dplyr)
library(multcompView)
library(tidyr)
library(purrr)
#---------------------------------------------------------
# 1. Function to run ANOVA + Tukey + extract letters
#---------------------------------------------------------
get_letters <- function(df) {
mod <- aov(ICR_Avail_P_avg ~ CA_typology, data = df)
tuk <- TukeyHSD(mod)
lets <- multcompLetters(tuk$CA_typology[, "p adj"])$Letters
data.frame(
CA_typology = names(lets),
Letters = lets,
stringsAsFactors = FALSE
)
}
#---------------------------------------------------------
# 2. Run statistics PER ENVIRONMENT
#---------------------------------------------------------
stats_df <- drivers_data %>%
filter(!is.na(CA_typology), !is.na(ICR_Avail_P_avg)) %>%
group_by(IQF_environment) %>%
nest() %>%
mutate(
letters = map(data, get_letters),
means = map(data, ~ .x %>%
group_by(CA_typology) %>%
summarise(mean_val = mean(ICR_Avail_P_avg, na.rm = TRUE),
.groups = "drop"))
) %>%
unnest(letters) %>%
left_join(
unnest(dplyr::select(., IQF_environment, means)),
by = c("IQF_environment", "CA_typology")
) %>%
mutate(mean_label = paste0("mean = ", round(mean_val, 2)))
#---------------------------------------------------------
# 3. Plot (letters + means fixed at Y)
#---------------------------------------------------------
ggplot(drivers_data, aes(x = CA_typology, y = ICR_Avail_P_avg)) +
geom_boxplot(outlier.shape = NA, fill = "lightgray") +
geom_jitter(width = 0.2, alpha = 0.6, color = "steelblue", size = 1.4) +
# Letters
geom_text(
data = stats_df,
aes(x = CA_typology, y = 110, label = Letters),
inherit.aes = FALSE,
fontface = "bold",
size = 3,
color = "darkred"
) +
# Means
geom_text(
data = stats_df,
aes(x = CA_typology, y = 100, label = mean_label),
inherit.aes = FALSE,
size = 3
) +
facet_wrap(~ IQF_environment) +
labs(
title = "ICR_Avail_P_avg by CA_typology (per environment)",
subtitle = "Letters = Tukey HSD within each environment",
x = "CA Typology",
y = "Available Phosphorus (ppm)"
) +
theme_minimal()
###Correlation: ICR_Org_C_avg vs. ICR_Avail_P_avg
First we check the values. I suspect they are toon high
library(dplyr)
library(ggplot2)
drivers_data_clean <- drivers_data %>%
filter(!is.na(ICR_Avail_P_avg))
ggplot(drivers_data_clean,
aes(x = ICR_Avail_P_avg)) +
stat_ecdf(geom = "step", size = 1.2, color = "blue") +
scale_x_continuous(
breaks = seq(
floor(min(drivers_data_clean$ICR_Avail_P_avg, na.rm = TRUE)),
ceiling(max(drivers_data_clean$ICR_Avail_P_avg, na.rm = TRUE)),
by = 15
)
) +
labs(
title = "Cumulative Distribution of Available P",
x = "ICR_Avail_P_avg",
y = "Cumulative Proportion"
) +
theme_minimal()
they still seem high, lets see the relation with the evironmetal index
calculated before. This is an index of the yield of all the seasons in
the block as compared to plots in the same seasson-district
library(dplyr)
library(ggplot2)
library(segmented)
# ===============================
# 1) Average env_index per block across seasons
# ===============================
env_block_avg <- CONV_CA_Ratio %>%
filter(!is.na(IQR_block), !is.na(env_index)) %>%
group_by(IQR_block) %>%
summarise(
env_index_avg = mean(env_index, na.rm = TRUE),
.groups = "drop"
)
# ===============================
# 2) Ensure drivers_data is block-level for P
# ===============================
p_block <- drivers_data %>%
filter(!is.na(IQR_block), !is.na(ICR_Avail_P_avg)) %>%
group_by(IQR_block) %>%
summarise(
ICR_Avail_P_avg = mean(ICR_Avail_P_avg, na.rm = TRUE),
.groups = "drop"
)
# ===============================
# 3) Join and restrict to P ≤ 70
# ===============================
plot_data <- env_block_avg %>%
left_join(p_block, by = "IQR_block") %>%
filter(!is.na(ICR_Avail_P_avg),
ICR_Avail_P_avg <= 70)
cat("Blocks used:", nrow(plot_data), "\n")## Blocks used: 442# ===============================
# 4) Fit linear model
# ===============================
lm_model <- lm(env_index_avg ~ ICR_Avail_P_avg,
data = plot_data)
# ===============================
# 5) Fit segmented (breakpoint) model
# ===============================
seg_model <- segmented(
lm_model,
seg.Z = ~ ICR_Avail_P_avg,
psi = list(ICR_Avail_P_avg = 30) # initial guess
)
summary(seg_model)##
## ***Regression Model with Segmented Relationship(s)***
##
## Call:
## segmented.lm(obj = lm_model, seg.Z = ~ICR_Avail_P_avg, psi = list(ICR_Avail_P_avg = 30))
##
## Estimated Break-Point(s):
## Est. St.Err
## psi1.ICR_Avail_P_avg 4.829 0.999
##
## Coefficients of the linear terms:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 0.48109 0.32650 1.473 0.141
## ICR_Avail_P_avg 0.10610 0.08338 1.272 0.204
## U1.ICR_Avail_P_avg -0.10491 0.08339 -1.258 NA
##
## Residual standard error: 0.2971 on 438 degrees of freedom
## Multiple R-Squared: 0.01508, Adjusted R-squared: 0.008334
##
## Boot restarting based on 8 samples. Last fit:
## Convergence attained in 1 iterations (rel. change 0)# Extract breakpoint
breakpoint <- seg_model$psi[2]
cat("Estimated Breakpoint:", round(breakpoint, 2), "\n")## Estimated Breakpoint: 4.83# Add predictions
plot_data$predicted <- predict(seg_model)
# ===============================
# 6) Plot
# ===============================
ggplot(plot_data,
aes(x = ICR_Avail_P_avg,
y = env_index_avg)) +
geom_point(alpha = 0.7,
color = "steelblue",
size = 2) +
geom_line(aes(y = predicted),
color = "red",
size = 1.2) +
geom_vline(xintercept = breakpoint,
linetype = "dashed",
color = "darkgreen",
size = 1) +
scale_x_continuous(breaks = seq(0, 70, by = 10)) +
labs(
title = "Breakpoint Model: Environmental Index vs Available P",
subtitle = paste0("Estimated breakpoint ≈ ",
round(breakpoint, 2)),
x = "ICR_Avail_P_avg",
y = "Mean env_index across seasons"
) +
theme_minimal(base_size = 14)
library(ggplot2)
library(dplyr)
library(purrr)
library(tidyr)
#---------------------------------------------------------
# 1. Fit linear models per environment
#---------------------------------------------------------
correlation_stats <- drivers_data %>%
filter(!is.na(ICR_Org_C_avg), !is.na(ICR_Avail_P_avg), !is.na(IQF_environment)) %>%
group_by(IQF_environment) %>%
nest() %>%
mutate(
model = map(data, ~ lm(ICR_Avail_P_avg ~ ICR_Org_C_avg, data = .x)),
p_value = map_dbl(model, ~ summary(.x)$coefficients[2, 4]),
r_squared = map_dbl(model, ~ summary(.x)$r.squared),
stat_label = paste0("p = ", signif(p_value, 3),
"\nR", "\u00B2", " = ", round(r_squared, 2))
)
#---------------------------------------------------------
# 2. Set annotation position
#---------------------------------------------------------
label_y <- max(drivers_data$ICR_Avail_P_avg, na.rm = TRUE) * 0.95
label_x <- max(drivers_data$ICR_Org_C_avg, na.rm = TRUE) * 0.95
correlation_stats <- correlation_stats %>%
mutate(x = label_x, y = label_y)
#---------------------------------------------------------
# 3. Plot: scatter + regression + p + R² + facet
#---------------------------------------------------------
ggplot(drivers_data, aes(x = ICR_Org_C_avg, y = ICR_Avail_P_avg)) +
geom_point(alpha = 0.6, color = "steelblue", size = 1.5) +
geom_smooth(method = "lm", se = TRUE, color = "darkred", linetype = "dashed") +
geom_text(
data = correlation_stats,
aes(x = x, y = y, label = stat_label),
inherit.aes = FALSE,
hjust = 1, vjust = 1,
size = 4
) +
facet_wrap(~ IQF_environment) +
labs(
title = "Correlation between Organic Carbon and Available Phosphorus",
subtitle = "Linear regression per environment with p-values and R\u00B2",
x = "ICR_Org_C_avg (Organic Carbon %)",
y = "ICR_Avail_P_avg (Phosphorus ppm)"
) +
theme_minimal()
CA Performance as affected by ICR_arable_land_avg
NOtes: no correlation when pooled together but a trant in 2 out of 4 environments. Moreover it is a hosehold variable that is ususally used to explain perfomrance of practices Conclusions: Keep it for drivers analysis
with avg_Y_ratio
library(ggplot2)
# 1. Fit linear regression model
model <- lm(avg_Y_ratio ~ ICR_arable_land_avg, data = drivers_data)
# 2. Extract p-value from model summary
p_val <- summary(model)$coefficients[2, 4] # p-value for slope
# 3. Create plot
ggplot(drivers_data, aes(x = ICR_arable_land_avg, y = avg_Y_ratio)) +
geom_point(alpha = 0.6, color = "steelblue") +
geom_smooth(method = "lm", se = TRUE, color = "darkred", linetype = "dashed") +
labs(
title = "ICR_arable_land_avg vs. avg_Y_ratio",
x = "Arable Land (%)",
y = "Average Yield Ratio"
) +
annotate("text",
x = Inf, y = Inf, hjust = 1.1, vjust = 1.5,
label = paste0("p = ", signif(p_val, 3)),
size = 5, color = "black") +
theme_minimal()
### In facet by environment
library(ggplot2)
library(dplyr)
library(purrr)
library(tidyr)
#----------------------------------------------
# 1. Fit separate linear models per environment
#----------------------------------------------
stats_df <- drivers_data %>%
filter(!is.na(ICR_arable_land_avg), !is.na(avg_Y_ratio), !is.na(IQF_environment)) %>%
group_by(IQF_environment) %>%
nest() %>%
mutate(
model = map(data, ~ lm(avg_Y_ratio ~ ICR_arable_land_avg, data = .x)),
p_value = map_dbl(model, ~ summary(.x)$coefficients[2, 4]),
p_label = paste0("p = ", signif(p_value, 3))
) %>%
dplyr::select(IQF_environment, p_label)
#----------------------------------------------
# 2. Prepare annotation positions
#----------------------------------------------
label_y <- max(drivers_data$avg_Y_ratio, na.rm = TRUE) * 0.95
label_x <- max(drivers_data$ICR_arable_land_avg, na.rm = TRUE) * 0.95
stats_df <- stats_df %>%
mutate(x = label_x, y = label_y)
#----------------------------------------------
# 3. Plot
#----------------------------------------------
ggplot(drivers_data, aes(x = ICR_arable_land_avg, y = avg_Y_ratio)) +
geom_point(alpha = 0.6, color = "steelblue") +
geom_smooth(method = "lm", se = TRUE, color = "darkred", linetype = "dashed") +
geom_text(
data = stats_df,
aes(x = x, y = y, label = p_label),
inherit.aes = FALSE,
size = 5,
hjust = 1,
vjust = 1
) +
facet_wrap(~ IQF_environment) +
labs(
title = "ICR_arable_land_avg vs. avg_Y_ratio",
subtitle = "Linear regressions and p-values shown per environment",
x = "Arable land (%)",
y = "Average Yield Ratio"
) +
theme_minimal()
###ICR_arable_land_avg by CA_typology
library(ggplot2)
library(dplyr)
library(multcompView)
library(tidyr)
library(purrr)
get_letters <- function(df) {
mod <- aov(ICR_arable_land_avg ~ CA_typology, data = df)
tuk <- TukeyHSD(mod)
lets <- multcompLetters(tuk$CA_typology[, "p adj"])$Letters
data.frame(
CA_typology = names(lets),
Letters = lets,
stringsAsFactors = FALSE
)
}
stats_df <- drivers_data %>%
filter(!is.na(CA_typology), !is.na(ICR_arable_land_avg)) %>%
group_by(IQF_environment) %>%
nest() %>%
mutate(
letters = map(data, get_letters),
means = map(data, ~ .x %>%
group_by(CA_typology) %>%
summarise(mean_val = mean(ICR_arable_land_avg, na.rm = TRUE),
.groups = "drop"))
) %>%
unnest(letters) %>%
left_join(
unnest(dplyr::select(., IQF_environment, means)),
by = c("IQF_environment", "CA_typology")
) %>%
mutate(mean_label = paste0("mean = ", round(mean_val, 1)))
ggplot(drivers_data, aes(x = CA_typology, y = ICR_arable_land_avg)) +
geom_boxplot(outlier.shape = NA, fill = "lightgray") +
geom_jitter(width = 0.2, alpha = 0.6, color = "steelblue", size = 1.4) +
geom_text(
data = stats_df,
aes(x = CA_typology, y = 90, label = Letters),
inherit.aes = FALSE,
fontface = "bold",
size = 3,
color = "darkred"
) +
geom_text(
data = stats_df,
aes(x = CA_typology, y = 82, label = mean_label),
inherit.aes = FALSE,
size = 3
) +
facet_wrap(~ IQF_environment) +
labs(
title = "ICR_arable_land_avg by CA_typology (per environment)",
subtitle = "Letters = Tukey HSD within each environment",
x = "CA Typology",
y = "Arable land (%)"
) +
theme_minimal()
###Correlation: ICR_Org_C_avg vs ICR_arable_land_avg
library(ggplot2)
library(dplyr)
library(purrr)
library(tidyr)
correlation_stats <- drivers_data %>%
filter(!is.na(ICR_Org_C_avg), !is.na(ICR_arable_land_avg), !is.na(IQF_environment)) %>%
group_by(IQF_environment) %>%
nest() %>%
mutate(
model = map(data, ~ lm(ICR_arable_land_avg ~ ICR_Org_C_avg, data = .x)),
p_value = map_dbl(model, ~ summary(.x)$coefficients[2, 4]),
r_squared = map_dbl(model, ~ summary(.x)$r.squared),
stat_label = paste0("p = ", signif(p_value, 3),
"\nR", "\u00B2", " = ", round(r_squared, 2))
)
label_y <- max(drivers_data$ICR_arable_land_avg, na.rm = TRUE) * 0.95
label_x <- max(drivers_data$ICR_Org_C_avg, na.rm = TRUE) * 0.95
correlation_stats <- correlation_stats %>%
mutate(x = label_x, y = label_y)
ggplot(drivers_data, aes(x = ICR_Org_C_avg, y = ICR_arable_land_avg)) +
geom_point(alpha = 0.6, color = "steelblue", size = 1.5) +
geom_smooth(method = "lm", se = TRUE, color = "darkred", linetype = "dashed") +
geom_text(
data = correlation_stats,
aes(x = x, y = y, label = stat_label),
inherit.aes = FALSE,
hjust = 1, vjust = 1,
size = 4
) +
facet_wrap(~ IQF_environment) +
labs(
title = "Correlation between Organic Carbon and Arable Land",
subtitle = "Linear regression per environment with p-values and R\u00B2",
x = "ICR_Org_C_avg (Organic Carbon %)",
y = "ICR_arable_land_avg (Arable land %)"
) +
theme_minimal()
## CA Performance as affected by ICF_Alt_avg
NOtes: there seems to be no correlation when pooled together and even less when split by environment. Environment clasiffication is mainly based on altitude so if we are using em=nvironments as driver, altitude can be eliminated Conclusions: Do not use further in drivers analysis
###All together
library(ggplot2)
# 1. Fit linear regression model
model <- lm(avg_Y_ratio ~ ICF_Alt_avg, data = drivers_data)
# 2. Extract p-value from model summary
p_val <- summary(model)$coefficients[2, 4] # p-value for slope
# 3. Create plot
ggplot(drivers_data, aes(x = ICF_Alt_avg, y = avg_Y_ratio)) +
geom_point(alpha = 0.6, color = "steelblue") +
geom_smooth(method = "lm", se = TRUE, color = "darkred", linetype = "dashed") +
labs(
title = "ICF_Alt_avg vs. avg_Y_ratio",
x = "Altitude (%)",
y = "Average Yield Ratio"
) +
annotate("text",
x = Inf, y = Inf, hjust = 1.1, vjust = 1.5,
label = paste0("p = ", signif(p_val, 3)),
size = 5, color = "black") +
theme_minimal()
with avg_Y_ratio
library(ggplot2)
library(dplyr)
library(purrr)
library(tidyr)
#----------------------------------------------
# 1. Fit separate linear models per environment
#----------------------------------------------
stats_df <- drivers_data %>%
filter(!is.na(ICF_Alt_avg), !is.na(avg_Y_ratio), !is.na(IQF_environment)) %>%
group_by(IQF_environment) %>%
nest() %>%
mutate(
model = map(data, ~ lm(avg_Y_ratio ~ ICF_Alt_avg, data = .x)),
p_value = map_dbl(model, ~ summary(.x)$coefficients[2, 4]),
p_label = paste0("p = ", signif(p_value, 3))
) %>%
dplyr::select(IQF_environment, p_label)
#----------------------------------------------
# 2. Annotation positions
#----------------------------------------------
label_y <- max(drivers_data$avg_Y_ratio, na.rm = TRUE) * 0.95
label_x <- max(drivers_data$ICF_Alt_avg, na.rm = TRUE) * 0.95
stats_df <- stats_df %>%
mutate(x = label_x, y = label_y)
#----------------------------------------------
# 3. Plot
#----------------------------------------------
ggplot(drivers_data, aes(x = ICF_Alt_avg, y = avg_Y_ratio)) +
geom_point(alpha = 0.6, color = "steelblue") +
geom_smooth(method = "lm", se = TRUE, color = "darkred", linetype = "dashed") +
geom_text(
data = stats_df,
aes(x = x, y = y, label = p_label),
inherit.aes = FALSE,
size = 5,
hjust = 1,
vjust = 1
) +
facet_wrap(~ IQF_environment) +
labs(
title = "ICF_Alt_avg vs. avg_Y_ratio",
subtitle = "Linear regressions and p-values per environment",
x = "Altitude (m)",
y = "Average Yield Ratio"
) +
theme_minimal()
###ICF_Alt_avg by CA_typology
library(ggplot2)
library(dplyr)
library(multcompView)
library(tidyr)
library(purrr)
# Function for letters
get_letters <- function(df) {
mod <- aov(ICF_Alt_avg ~ CA_typology, data = df)
tuk <- TukeyHSD(mod)
lets <- multcompLetters(tuk$CA_typology[, "p adj"])$Letters
data.frame(
CA_typology = names(lets),
Letters = lets,
stringsAsFactors = FALSE
)
}
# Run stats by environment
stats_df <- drivers_data %>%
filter(!is.na(CA_typology), !is.na(ICF_Alt_avg)) %>%
group_by(IQF_environment) %>%
nest() %>%
mutate(
letters = map(data, get_letters),
means = map(data, ~ .x %>%
group_by(CA_typology) %>%
summarise(mean_val = mean(ICF_Alt_avg, na.rm = TRUE),
.groups = "drop"))
) %>%
unnest(letters) %>%
left_join(
unnest(dplyr::select(., IQF_environment, means)),
by = c("IQF_environment", "CA_typology")
) %>%
mutate(mean_label = paste0("mean = ", round(mean_val, 0)))
# Plot
ggplot(drivers_data, aes(x = CA_typology, y = ICF_Alt_avg)) +
geom_boxplot(outlier.shape = NA, fill = "lightgray") +
geom_jitter(width = 0.2, alpha = 0.6, color = "steelblue", size = 1.4) +
geom_text(
data = stats_df,
aes(x = CA_typology, y = 2000, label = Letters),
inherit.aes = FALSE,
fontface = "bold",
size = 3,
color = "darkred"
) +
geom_text(
data = stats_df,
aes(x = CA_typology, y = 1800, label = mean_label),
inherit.aes = FALSE,
size = 3
) +
facet_wrap(~ IQF_environment) +
labs(
title = "ICF_Alt_avg by CA_typology (per environment)",
subtitle = "Tukey HSD: significant differences by letter",
x = "CA Typology",
y = "Altitude (m)"
) +
theme_minimal()
CA Performance as affected by ICR_rainfall_avg
NOtes: Rain seems to partially explain CA performance, with a cuadratice-model response when pooled all environmetns together, with 4 mm per day as optimumn. However, when split by environments, linear models can fit as good Conclusions: Keep, consider cuadratic relation and not linear if not grouped within environments ###Simple regression: ICR_rainfall_avg vs avg_Y_ratio
library(ggplot2)
# 1. Fit linear regression model
model <- lm(avg_Y_ratio ~ ICR_rainfall_avg, data = drivers_data)
# 2. Extract p-value and R-squared from model summary
p_val <- summary(model)$coefficients[2, 4]
r2 <- summary(model)$r.squared
# 3. Create plot
ggplot(drivers_data, aes(x = ICR_rainfall_avg, y = avg_Y_ratio)) +
geom_point(alpha = 0.6, color = "steelblue") +
geom_smooth(method = "lm", se = TRUE, color = "darkred", linetype = "dashed") +
labs(
title = "ICR_rainfall_avg vs. avg_Y_ratio",
x = "Rainfall (avg, mm)",
y = "Average Yield Ratio"
) +
annotate("text",
x = Inf, y = Inf, hjust = 1.1, vjust = 1.5,
label = paste0("p = ", signif(p_val, 3),
"\nR^2 = ", round(r2, 2)),
size = 5, color = "black") +
theme_minimal()
### What wbout tieh a cuadratic model?
library(ggplot2)
# 1. Fit quadratic model
model <- lm(avg_Y_ratio ~ ICR_rainfall_avg + I(ICR_rainfall_avg^2),
data = drivers_data)
# 2. Extract p-value and R²
p_val <- summary(model)$coefficients["I(ICR_rainfall_avg^2)", 4]
r2 <- summary(model)$r.squared
# 3. Create plot with p and R²
ggplot(drivers_data, aes(x = ICR_rainfall_avg, y = avg_Y_ratio)) +
geom_point(alpha = 0.6, color = "steelblue") +
geom_smooth(
method = "lm",
formula = y ~ x + I(x^2),
se = TRUE,
color = "darkred",
linetype = "dashed"
) +
labs(
title = "ICR_rainfall_avg vs. avg_Y_ratio (Quadratic)",
x = "Rainfall (avg, mm)",
y = "Average Yield Ratio"
) +
annotate("text",
x = Inf, y = Inf, hjust = 1.1, vjust = 1.5,
label = paste0("p = ", signif(p_val, 3),
"\nR^2 = ", round(r2, 2)),
size = 5)
###Faceted regression: ICR_rainfall_avg vs avg_Y_ratio by environment
library(ggplot2)
library(dplyr)
library(purrr)
library(tidyr)
# 1. Fit separate linear models per environment
stats_df <- drivers_data %>%
filter(!is.na(ICR_rainfall_avg), !is.na(avg_Y_ratio), !is.na(IQF_environment)) %>%
group_by(IQF_environment) %>%
nest() %>%
mutate(
model = map(data, ~ lm(avg_Y_ratio ~ ICR_rainfall_avg, data = .x)),
p_value = map_dbl(model, ~ summary(.x)$coefficients[2, 4]),
p_label = paste0("p = ", signif(p_value, 3))
) %>%
dplyr::select(IQF_environment, p_label)
# 2. Annotation positions
label_y <- max(drivers_data$avg_Y_ratio, na.rm = TRUE) * 0.95
label_x <- max(drivers_data$ICR_rainfall_avg, na.rm = TRUE) * 0.95
stats_df <- stats_df %>%
mutate(x = label_x, y = label_y)
# 3. Plot
ggplot(drivers_data, aes(x = ICR_rainfall_avg, y = avg_Y_ratio)) +
geom_point(alpha = 0.6, color = "steelblue") +
geom_smooth(method = "lm", se = TRUE, color = "darkred", linetype = "dashed") +
geom_text(
data = stats_df,
aes(x = x, y = y, label = p_label),
inherit.aes = FALSE,
size = 5,
hjust = 1,
vjust = 1
) +
facet_wrap(~ IQF_environment) +
labs(
title = "ICR_rainfall_avg vs. avg_Y_ratio",
subtitle = "Linear regressions and p-values per environment",
x = "Rainfall (avg, mm)",
y = "Average Yield Ratio"
) +
theme_minimal()
### same as avobe but with cuadratic model
library(ggplot2)
library(dplyr)
library(purrr)
library(tidyr)
# 1. Fit quadratic models per environment
stats_df <- drivers_data %>%
filter(!is.na(ICR_rainfall_avg), !is.na(avg_Y_ratio), !is.na(IQF_environment)) %>%
group_by(IQF_environment) %>%
nest() %>%
mutate(
model = map(data, ~ lm(avg_Y_ratio ~ ICR_rainfall_avg + I(ICR_rainfall_avg^2), data = .x)),
p_value = map_dbl(model, ~ summary(.x)$coefficients["I(ICR_rainfall_avg^2)", 4]),
r_squared = map_dbl(model, ~ summary(.x)$r.squared),
# IMPORTANT: use R^2 (not R²)
stat_label = paste0(
"p = ", signif(p_value, 3),
"\nR^2 = ", round(r_squared, 2)
)
) %>%
dplyr::select(IQF_environment, stat_label)
stats_df <- stats_df %>%
mutate(x = label_x, y = label_y)
# 3. Plot
ggplot(drivers_data, aes(x = ICR_rainfall_avg, y = avg_Y_ratio)) +
geom_point(alpha = 0.6, color = "steelblue") +
geom_smooth(
method = "lm",
formula = y ~ x + I(x^2),
se = TRUE,
color = "darkred",
linetype = "dashed"
) +
geom_text(
data = stats_df,
aes(x = x, y = y, label = stat_label),
inherit.aes = FALSE,
size = 5,
hjust = 1,
vjust = 1
) +
facet_wrap(~ IQF_environment) +
labs(
title = "ICR_rainfall_avg vs. avg_Y_ratio (Quadratic)",
subtitle = "Quadratic regressions per environment with p-values and R^2",
x = "Rainfall (avg, mm)",
y = "Average Yield Ratio"
)
###Boxplot: ICR_rainfall_avg by CA_typology (with Tukey HSD letters)
library(ggplot2)
library(dplyr)
library(multcompView)
library(tidyr)
library(purrr)
# Function to extract Tukey letters
get_letters <- function(df) {
mod <- aov(ICR_rainfall_avg ~ CA_typology, data = df)
tuk <- TukeyHSD(mod)
lets <- multcompLetters(tuk$CA_typology[, "p adj"])$Letters
data.frame(
CA_typology = names(lets),
Letters = lets,
stringsAsFactors = FALSE
)
}
# Run stats by environment
stats_df <- drivers_data %>%
filter(!is.na(CA_typology), !is.na(ICR_rainfall_avg)) %>%
group_by(IQF_environment) %>%
nest() %>%
mutate(
letters = map(data, get_letters),
means = map(data, ~ .x %>%
group_by(CA_typology) %>%
summarise(mean_val = mean(ICR_rainfall_avg, na.rm = TRUE),
.groups = "drop"))
) %>%
unnest(letters) %>%
left_join(
unnest(dplyr::select(., IQF_environment, means)),
by = c("IQF_environment", "CA_typology")
) %>%
mutate(mean_label = paste0("mean = ", round(mean_val, 2)))
# Plot
ggplot(drivers_data, aes(x = CA_typology, y = ICR_rainfall_avg)) +
geom_boxplot(outlier.shape = NA, fill = "lightgray") +
geom_jitter(width = 0.2, alpha = 0.6, color = "steelblue", size = 1.4) +
geom_text(
data = stats_df,
aes(x = CA_typology, y = 7, label = Letters),
inherit.aes = FALSE,
fontface = "bold",
size = 3,
color = "darkred"
) +
geom_text(
data = stats_df,
aes(x = CA_typology, y = 6.2, label = mean_label),
inherit.aes = FALSE,
size = 3
) +
facet_wrap(~ IQF_environment) +
labs(
title = "ICR_rainfall_avg by CA_typology (per environment)",
subtitle = "Tukey HSD: significant differences by letter",
x = "CA Typology",
y = "Rainfall (avg, mm)"
) +
theme_minimal()
CA Performance as affected by Comp_amount_corr_quali
NOtes: It seems to have a significant relation with both Y_relative and typology. Wether pooled all together of per environment, but when split by environment, all show the sam esign of trend, but is clear and significant only for mid-dry and transition Conclusions: Keep for the drivers analysis ###Simple regression: Comp_amount_corr_quali vs avg_Y_ratio
library(ggplot2)
# Filter data
df_filt <- drivers_data %>%
dplyr::filter(!is.na(Comp_amount_corr_quali),
!is.na(avg_Y_ratio),
Comp_amount_corr_quali <= 15)
# 1. Fit linear regression model
model <- lm(avg_Y_ratio ~ Comp_amount_corr_quali, data = df_filt)
# 2. Extract p-value and R^2
p_val <- summary(model)$coefficients[2, 4]
r2 <- summary(model)$r.squared
# 3. Create plot
ggplot(df_filt, aes(x = Comp_amount_corr_quali, y = avg_Y_ratio)) +
geom_point(alpha = 0.6, color = "steelblue") +
geom_smooth(method = "lm", se = TRUE, color = "darkred", linetype = "dashed") +
labs(
title = "Comp_amount_corr_quali vs. avg_Y_ratio",
x = "Compost amount (corrected, ≤15)",
y = "Average Yield Ratio"
) +
annotate("text",
x = Inf, y = Inf, hjust = 1.1, vjust = 1.5,
label = paste0("p = ", signif(p_val, 3),
"\nR^2 = ", round(r2, 2)),
size = 5) +
theme_minimal()
###Linear regression per environment (facet + R² + p, compost ≤ 15)
library(ggplot2)
library(dplyr)
library(purrr)
library(tidyr)
df_filt <- drivers_data %>%
dplyr::filter(!is.na(Comp_amount_corr_quali),
!is.na(avg_Y_ratio),
!is.na(IQF_environment),
Comp_amount_corr_quali <= 15)
#----------------------------------------------
# 1. Fit models per environment
#----------------------------------------------
stats_df <- df_filt %>%
group_by(IQF_environment) %>%
nest() %>%
mutate(
model = map(data, ~ lm(avg_Y_ratio ~ Comp_amount_corr_quali, data = .x)),
p_value = map_dbl(model, ~ summary(.x)$coefficients[2, 4]),
r_squared = map_dbl(model, ~ summary(.x)$r.squared),
stat_label = paste0("p = ", signif(p_value, 3),
"\nR^2 = ", round(r_squared, 2))
) %>%
dplyr::select(IQF_environment, stat_label)
#----------------------------------------------
# 2. Annotation positions
#----------------------------------------------
label_y <- max(df_filt$avg_Y_ratio, na.rm = TRUE) * 0.95
label_x <- max(df_filt$Comp_amount_corr_quali, na.rm = TRUE) * 0.95
stats_df <- stats_df %>%
mutate(x = label_x, y = label_y)
#----------------------------------------------
# 3. Plot
#----------------------------------------------
ggplot(df_filt, aes(x = Comp_amount_corr_quali, y = avg_Y_ratio)) +
geom_point(alpha = 0.6, color = "steelblue") +
geom_smooth(method = "lm", se = TRUE, color = "darkred", linetype = "dashed") +
geom_text(
data = stats_df,
aes(x = x, y = y, label = stat_label),
inherit.aes = FALSE,
size = 5,
hjust = 1,
vjust = 1
) +
facet_wrap(~ IQF_environment) +
labs(
title = "Comp_amount_corr_quali vs. avg_Y_ratio",
subtitle = "Linear regressions per environment (compost ≤ 15)",
x = "Compost amount (corrected)",
y = "Average Yield Ratio"
) +
theme_minimal()
###Boxplot by CA_typology (Tukey HSD, compost ≤ 15)
library(ggplot2)
library(dplyr)
library(multcompView)
library(tidyr)
library(purrr)
df_filt <- drivers_data %>%
dplyr::filter(!is.na(CA_typology),
!is.na(Comp_amount_corr_quali),
!is.na(IQF_environment),
Comp_amount_corr_quali <= 15)
get_letters <- function(df) {
mod <- aov(Comp_amount_corr_quali ~ CA_typology, data = df)
tuk <- TukeyHSD(mod)
lets <- multcompLetters(tuk$CA_typology[, "p adj"])$Letters
data.frame(
CA_typology = names(lets),
Letters = lets,
stringsAsFactors = FALSE
)
}
stats_df <- df_filt %>%
group_by(IQF_environment) %>%
nest() %>%
mutate(
letters = map(data, get_letters),
means = map(data, ~ .x %>%
group_by(CA_typology) %>%
summarise(mean_val = mean(Comp_amount_corr_quali, na.rm = TRUE),
.groups = "drop"))
) %>%
unnest(letters) %>%
left_join(
unnest(dplyr::select(., IQF_environment, means)),
by = c("IQF_environment", "CA_typology")
) %>%
mutate(mean_label = paste0("mean = ", round(mean_val, 2)))
ggplot(df_filt, aes(x = CA_typology, y = Comp_amount_corr_quali)) +
geom_boxplot(outlier.shape = NA, fill = "lightgray") +
geom_jitter(width = 0.2, alpha = 0.6, color = "steelblue", size = 1.4) +
geom_text(
data = stats_df,
aes(x = CA_typology,
y = max(df_filt$Comp_amount_corr_quali, na.rm = TRUE) * 0.95,
label = Letters),
inherit.aes = FALSE,
fontface = "bold",
size = 3,
color = "darkred"
) +
geom_text(
data = stats_df,
aes(x = CA_typology,
y = max(df_filt$Comp_amount_corr_quali, na.rm = TRUE) * 0.85,
label = mean_label),
inherit.aes = FALSE,
size = 3
) +
facet_wrap(~ IQF_environment) +
labs(
title = "Comp_amount_corr_quali by CA_typology (per environment)",
subtitle = "Tukey HSD (compost ≤ 15)",
x = "CA Typology",
y = "Compost amount (corrected)"
) +
theme_minimal()
## CA Performance as affected by ICF_planting_date
NOtes: Conclusions: ###Simple regression: ICR_rainfall_avg ICF_planting_date
library(ggplot2)
library(dplyr)
# Filter data
df_filt <- drivers_data %>%
dplyr::filter(!is.na(ICF_planting_date),
!is.na(avg_Y_ratio))
# 1. Fit linear regression model
model <- lm(avg_Y_ratio ~ ICF_planting_date, data = df_filt)
# 2. Extract p-value and R^2
p_val <- summary(model)$coefficients[2, 4]
r2 <- summary(model)$r.squared
# 3. Create plot
ggplot(df_filt, aes(x = ICF_planting_date, y = avg_Y_ratio)) +
geom_point(alpha = 0.6, color = "steelblue") +
geom_smooth(method = "lm", se = TRUE, color = "darkred", linetype = "dashed") +
labs(
title = "ICF_planting_date vs. avg_Y_ratio",
x = "Planting date",
y = "Average Yield Ratio"
) +
annotate("text",
x = Inf, y = Inf, hjust = 1.1, vjust = 1.5,
label = paste0("p = ", signif(p_val, 3),
"\nR^2 = ", round(r2, 2)),
size = 5) +
theme_minimal()
###Linear regression per environment (facet + R^2 + p)
library(ggplot2)
library(dplyr)
library(purrr)
library(tidyr)
df_filt <- drivers_data %>%
dplyr::filter(!is.na(ICF_planting_date),
!is.na(avg_Y_ratio),
!is.na(IQF_environment))
#----------------------------------------------
# 1. Fit models per environment
#----------------------------------------------
stats_df <- df_filt %>%
group_by(IQF_environment) %>%
nest() %>%
mutate(
model = map(data, ~ lm(avg_Y_ratio ~ ICF_planting_date, data = .x)),
p_value = map_dbl(model, ~ summary(.x)$coefficients[2, 4]),
r_squared = map_dbl(model, ~ summary(.x)$r.squared),
stat_label = paste0("p = ", signif(p_value, 3),
"\nR^2 = ", round(r_squared, 2))
) %>%
dplyr::select(IQF_environment, stat_label)
#----------------------------------------------
# 2. Annotation positions
#----------------------------------------------
label_y <- max(df_filt$avg_Y_ratio, na.rm = TRUE) * 0.95
label_x <- max(df_filt$ICF_planting_date, na.rm = TRUE) * 0.95
stats_df <- stats_df %>%
mutate(x = label_x, y = label_y)
#----------------------------------------------
# 3. Plot
#----------------------------------------------
ggplot(df_filt, aes(x = ICF_planting_date, y = avg_Y_ratio)) +
geom_point(alpha = 0.6, color = "steelblue") +
geom_smooth(method = "lm", se = TRUE, color = "darkred", linetype = "dashed") +
geom_text(
data = stats_df,
aes(x = x, y = y, label = stat_label),
inherit.aes = FALSE,
size = 5,
hjust = 1,
vjust = 1
) +
facet_wrap(~ IQF_environment) +
labs(
title = "ICF_planting_date vs. avg_Y_ratio",
subtitle = "Linear regressions per environment",
x = "Planting date",
y = "Average Yield Ratio"
) +
theme_minimal()
###Boxplot by CA_typology (Tukey HSD)
library(ggplot2)
library(dplyr)
library(multcompView)
library(tidyr)
library(purrr)
df_filt <- drivers_data %>%
dplyr::filter(!is.na(CA_typology),
!is.na(ICF_planting_date),
!is.na(IQF_environment))
get_letters <- function(df) {
mod <- aov(ICF_planting_date ~ CA_typology, data = df)
tuk <- TukeyHSD(mod)
lets <- multcompLetters(tuk$CA_typology[, "p adj"])$Letters
data.frame(
CA_typology = names(lets),
Letters = lets,
stringsAsFactors = FALSE
)
}
stats_df <- df_filt %>%
group_by(IQF_environment) %>%
nest() %>%
mutate(
letters = map(data, get_letters),
means = map(data, ~ .x %>%
group_by(CA_typology) %>%
summarise(mean_val = mean(ICF_planting_date, na.rm = TRUE),
.groups = "drop"))
) %>%
unnest(letters) %>%
left_join(
unnest(dplyr::select(., IQF_environment, means)),
by = c("IQF_environment", "CA_typology")
) %>%
mutate(mean_label = paste0("mean = ", round(mean_val, 1)))
ggplot(df_filt, aes(x = CA_typology, y = ICF_planting_date)) +
geom_boxplot(outlier.shape = NA, fill = "lightgray") +
geom_jitter(width = 0.2, alpha = 0.6, color = "steelblue", size = 1.4) +
geom_text(
data = stats_df,
aes(x = CA_typology,
y = max(df_filt$ICF_planting_date, na.rm = TRUE) * 0.95,
label = Letters),
inherit.aes = FALSE,
fontface = "bold",
size = 3,
color = "darkred"
) +
geom_text(
data = stats_df,
aes(x = CA_typology,
y = max(df_filt$ICF_planting_date, na.rm = TRUE) * 0.85,
label = mean_label),
inherit.aes = FALSE,
size = 3
) +
facet_wrap(~ IQF_environment) +
labs(
title = "ICF_planting_date by CA_typology (per environment)",
subtitle = "Tukey HSD",
x = "CA Typology",
y = "Planting date"
) +
theme_minimal()
###How does delay in planting relate to rain over the season?
library(ggplot2)
library(dplyr)
library(tidyr)
library(purrr)
# 1. Filter valid data
df_filt <- drivers_data %>%
dplyr::filter(
!is.na(ICF_planting_date),
!is.na(ICR_rainfall_avg),
!is.na(IQF_environment)
)
# 2. Fit linear model per environment
stats_df <- df_filt %>%
group_by(IQF_environment) %>%
nest() %>%
mutate(
model = map(data, ~ lm(ICR_rainfall_avg ~ ICF_planting_date, data = .x)),
p_value = map_dbl(model, ~ summary(.x)$coefficients[2, 4]),
r_squared = map_dbl(model, ~ summary(.x)$r.squared),
stat_label = paste0("p = ", signif(p_value, 3),
"\nR^2 = ", round(r_squared, 2))
) %>%
dplyr::select(IQF_environment, stat_label)
# 3. Annotation position
label_y <- max(df_filt$ICR_rainfall_avg, na.rm = TRUE) * 0.95
label_x <- max(df_filt$ICF_planting_date, na.rm = TRUE) * 0.95
stats_df <- stats_df %>%
mutate(x = label_x, y = label_y)
# 4. Plot
ggplot(df_filt, aes(x = ICF_planting_date, y = ICR_rainfall_avg)) +
geom_point(alpha = 0.6, color = "steelblue", size = 1.5) +
geom_smooth(method = "lm", se = TRUE, color = "darkred", linetype = "dashed") +
geom_text(
data = stats_df,
aes(x = x, y = y, label = stat_label),
inherit.aes = FALSE,
size = 4,
hjust = 1,
vjust = 1
) +
facet_wrap(~ IQF_environment) +
labs(
title = "ICF_planting_date vs. ICR_rainfall_avg",
subtitle = "Linear regression per environment (p-value and R^2 shown)",
x = "Planting date",
y = "Rainfall average (mm)"
) +
theme_minimal()
CA Performance as affected by slope of the plot
library(ggplot2)
library(dplyr)
# Filter data
df_filt <- drivers_data %>%
filter(!is.na(Slope),
!is.na(avg_Y_ratio))
# 1. Fit linear regression model
model <- lm(avg_Y_ratio ~ Slope, data = df_filt)
# 2. Extract p-value and R^2
p_val <- summary(model)$coefficients[2, 4]
r2 <- summary(model)$r.squared
# 3. Create plot
ggplot(df_filt, aes(x = Slope, y = avg_Y_ratio)) +
geom_point(alpha = 0.6, color = "steelblue") +
geom_smooth(method = "lm", se = TRUE,
color = "darkred", linetype = "dashed") +
labs(
title = "Slope vs. avg_Y_ratio",
x = "Slope",
y = "Average Yield Ratio"
) +
annotate("text",
x = Inf, y = Inf,
hjust = 1.1, vjust = 1.5,
label = paste0("p = ", signif(p_val, 3),
"\nR^2 == ", round(r2, 2)),
parse = TRUE,
size = 5) +
theme_minimal()
###Linear regression per environment using Slope
library(ggplot2)
library(dplyr)
library(purrr)
library(tidyr)
df_filt <- drivers_data %>%
filter(!is.na(Slope),
!is.na(avg_Y_ratio),
!is.na(IQF_environment))
#----------------------------------------------
# 1. Fit models per environment
#----------------------------------------------
stats_df <- df_filt %>%
group_by(IQF_environment) %>%
nest() %>%
mutate(
model = map(data, ~ lm(avg_Y_ratio ~ Slope, data = .x)),
p_value = map_dbl(model, ~ summary(.x)$coefficients[2, 4]),
r_squared = map_dbl(model, ~ summary(.x)$r.squared),
stat_label = paste0("p = ", signif(p_value, 3),
"\nR^2 = ", round(r_squared, 2))
) %>%
dplyr::select(IQF_environment, stat_label)
#----------------------------------------------
# 2. Annotation positions
#----------------------------------------------
label_y <- max(df_filt$avg_Y_ratio, na.rm = TRUE) * 0.95
label_x <- max(df_filt$Slope, na.rm = TRUE) * 0.95
stats_df <- stats_df %>%
mutate(x = label_x, y = label_y)
#----------------------------------------------
# 3. Plot
#----------------------------------------------
ggplot(df_filt, aes(x = Slope, y = avg_Y_ratio)) +
geom_point(alpha = 0.6, color = "steelblue") +
geom_smooth(method = "lm", se = TRUE,
color = "darkred", linetype = "dashed") +
geom_text(
data = stats_df,
aes(x = x, y = y, label = stat_label),
inherit.aes = FALSE,
size = 5,
hjust = 1,
vjust = 1
) +
facet_wrap(~ IQF_environment) +
labs(
title = "Slope vs. avg_Y_ratio",
subtitle = "Linear regressions per environment",
x = "Slope",
y = "Average Yield Ratio"
) +
theme_minimal()
### Now slopes per CA_typology
library(ggplot2)
library(dplyr)
library(multcompView)
library(tidyr)
library(purrr)
df_filt <- drivers_data %>%
filter(!is.na(CA_typology),
!is.na(Slope),
!is.na(IQF_environment))
#----------------------------------------------
# Function to extract Tukey letters
#----------------------------------------------
get_letters <- function(df) {
mod <- aov(Slope ~ CA_typology, data = df)
tuk <- TukeyHSD(mod)
lets <- multcompLetters(tuk$CA_typology[, "p adj"])$Letters
data.frame(
CA_typology = names(lets),
Letters = lets,
stringsAsFactors = FALSE
)
}
#----------------------------------------------
# Compute letters and means
#----------------------------------------------
stats_df <- df_filt %>%
group_by(IQF_environment) %>%
nest() %>%
mutate(
letters = map(data, get_letters),
means = map(data, ~ .x %>%
group_by(CA_typology) %>%
summarise(mean_val = mean(Slope, na.rm = TRUE),
.groups = "drop"))
) %>%
unnest(letters) %>%
left_join(
unnest(dplyr::select(., IQF_environment, means)),
by = c("IQF_environment", "CA_typology")
) %>%
mutate(mean_label = paste0("mean = ", round(mean_val, 2)))
#----------------------------------------------
# Plot
#----------------------------------------------
ggplot(df_filt, aes(x = CA_typology, y = Slope)) +
geom_boxplot(outlier.shape = NA, fill = "lightgray") +
geom_jitter(width = 0.2, alpha = 0.6,
color = "steelblue", size = 1.4) +
geom_text(
data = stats_df,
aes(x = CA_typology,
y = max(df_filt$Slope, na.rm = TRUE) * 0.95,
label = Letters),
inherit.aes = FALSE,
fontface = "bold",
size = 3,
color = "darkred"
) +
geom_text(
data = stats_df,
aes(x = CA_typology,
y = max(df_filt$Slope, na.rm = TRUE) * 0.85,
label = mean_label),
inherit.aes = FALSE,
size = 3
) +
facet_wrap(~ IQF_environment) +
labs(
title = "Slope by CA_typology (per environment)",
subtitle = "ANOVA + Tukey HSD",
x = "CA Typology",
y = "Slope"
) +
theme_minimal()
##Factorial variables ## CA Performance as affected by IQR_TF_tubura_client NOtes: Conclusions: ###IQR_TF_tubura_client with Y_ratio
library(ggplot2)
library(dplyr)
# 1. Filter valid data
df_filt <- drivers_data %>%
dplyr::filter(
!is.na(avg_Y_ratio),
!is.na(IQR_TF_tubura_client),
!IQR_TF_tubura_client %in% c("no data", "---", "", "NA")
) %>%
dplyr::mutate(IQR_TF_tubura_client = droplevels(factor(IQR_TF_tubura_client))) # convert to factor, then drop unused levels
# 2. Run ANOVA
anova_model <- aov(avg_Y_ratio ~ IQR_TF_tubura_client, data = df_filt)
anova_p <- summary(anova_model)[[1]][["Pr(>F)"]][1]
# 3. Run Kruskal-Wallis test
kruskal_test <- kruskal.test(avg_Y_ratio ~ IQR_TF_tubura_client, data = df_filt)
kruskal_p <- kruskal_test$p.value
# 4. Build annotation label
test_label <- paste0(
"ANOVA p = ", signif(anova_p, 3),
"\nKruskal p = ", signif(kruskal_p, 3)
)
# 5. Calculate group means
mean_df <- df_filt %>%
group_by(IQR_TF_tubura_client) %>%
summarise(mean_val = mean(avg_Y_ratio, na.rm = TRUE), .groups = "drop") %>%
mutate(mean_label = paste0("mean = ", round(mean_val, 2)))
# 6. Plot boxplot with p-values and means
ggplot(df_filt, aes(x = IQR_TF_tubura_client, y = avg_Y_ratio)) +
geom_boxplot(outlier.shape = NA, fill = "lightgray") +
geom_jitter(width = 0.2, alpha = 0.6, color = "steelblue") +
geom_text(
data = mean_df,
aes(x = IQR_TF_tubura_client,
y = min(df_filt$avg_Y_ratio, na.rm = TRUE) * 0.95,
label = mean_label),
inherit.aes = FALSE,
size = 4
) +
annotate("text",
x = Inf, y = Inf,
hjust = 1.1, vjust = 1.5,
label = test_label,
size = 5) +
labs(
title = "avg_Y_ratio by IQR_TF_tubura_client",
subtitle = "ANOVA and Kruskal–Wallis test results shown",
x = "Tubura client status",
y = "Average Yield Ratio"
) +
theme_minimal()
###IQR_TF_tubura_client with CA_typology
library(ggplot2)
library(dplyr)
library(scales)
# Clean and prepare data
df_bar <- drivers_data %>%
filter(
!is.na(CA_typology),
!is.na(IQR_TF_tubura_client),
!CA_typology %in% c("no data", "---", "", "NA"),
!IQR_TF_tubura_client %in% c("no data", "---", "", "NA")
)
# Calculate N per category of explanatory variable
n_labels <- df_bar %>%
group_by(IQR_TF_tubura_client) %>%
summarise(n = n(), .groups = "drop")
# Plot with N annotations
ggplot(df_bar, aes(x = IQR_TF_tubura_client, fill = CA_typology)) +
geom_bar(position = "fill", color = "black") +
geom_text(
data = n_labels,
aes(x = IQR_TF_tubura_client, y = 1.02, label = paste0("N = ", n)),
inherit.aes = FALSE,
size = 4
) +
scale_y_continuous(labels = percent_format(), expand = expansion(mult = c(0, 0.1))) +
labs(
title = "CA_typology distribution by Tubura client status",
subtitle = "Proportions with total sample size per group",
x = "IQR_TF_tubura_client",
y = "Proportion of CA Typology",
fill = "CA Typology"
) +
theme_minimal()
# Build contingency table
tab <- table(df_bar$IQR_TF_tubura_client, df_bar$CA_typology)
# Run chi-squared test
chi_test <- chisq.test(tab)
# Print results
chi_test##
## Pearson's Chi-squared test
##
## data: tab
## X-squared = 1.443, df = 3, p-value = 0.6955# Fisher’s Exact Test (for categorical association)
fisher_test <- fisher.test(tab)
# Print results
fisher_test##
## Fisher's Exact Test for Count Data
##
## data: tab
## p-value = 0.673
## alternative hypothesis: two.sidedCA Performance as affected by IQR_plot_fert_quality
NOtes: Conclusions: ###IQR_plot_fert_quality with Y_ratio
library(ggplot2)
library(dplyr)
library(multcompView)
# 1. Filter and clean data
df_filt <- drivers_data %>%
filter(
!is.na(avg_Y_ratio),
!is.na(IQR_plot_fert_quality),
!IQR_plot_fert_quality %in% c("no data", "---", "", "NA")
) %>%
mutate(IQR_plot_fert_quality = factor(IQR_plot_fert_quality, levels = c("low_quality", "medium", "high")))
# 2. Run ANOVA
anova_model <- aov(avg_Y_ratio ~ IQR_plot_fert_quality, data = df_filt)
anova_p <- summary(anova_model)[[1]][["Pr(>F)"]][1]
# 3. Tukey HSD and Letters
tukey <- TukeyHSD(anova_model)
letters <- multcompLetters(tukey$`IQR_plot_fert_quality`[, "p adj"])$Letters
letters_df <- data.frame(
IQR_plot_fert_quality = names(letters),
Letters = letters,
stringsAsFactors = FALSE
)
# 4. Kruskal-Wallis test
kruskal_p <- kruskal.test(avg_Y_ratio ~ IQR_plot_fert_quality, data = df_filt)$p.value
test_label <- paste0(
"ANOVA p = ", signif(anova_p, 3),
"\nKruskal p = ", signif(kruskal_p, 3)
)
# 5. Group means
means_df <- df_filt %>%
group_by(IQR_plot_fert_quality) %>%
summarise(mean_val = mean(avg_Y_ratio, na.rm = TRUE), .groups = "drop") %>%
mutate(mean_label = paste0("mean = ", round(mean_val, 2)))
# 6. Combine means and letters
labels_df <- left_join(means_df, letters_df, by = "IQR_plot_fert_quality")
# 7. Plot
ggplot(df_filt, aes(x = IQR_plot_fert_quality, y = avg_Y_ratio)) +
geom_boxplot(outlier.shape = NA, fill = "lightgray") +
geom_jitter(width = 0.2, alpha = 0.6, color = "steelblue") +
# Add Tukey letters above boxes
geom_text(
data = labels_df,
aes(x = IQR_plot_fert_quality, y = max(df_filt$avg_Y_ratio, na.rm = TRUE) * 1.02, label = Letters),
inherit.aes = FALSE,
fontface = "bold", size = 5, color = "darkred"
) +
# Add means below boxes
geom_text(
data = labels_df,
aes(x = IQR_plot_fert_quality, y = min(df_filt$avg_Y_ratio, na.rm = TRUE) * 0.95, label = mean_label),
inherit.aes = FALSE,
size = 4
) +
# Add test label
annotate("text",
x = Inf, y = Inf,
hjust = 1.1, vjust = 1.5,
label = test_label,
size = 5) +
labs(
title = "avg_Y_ratio by IQR_plot_fert_quality",
subtitle = "Tukey HSD letters, means, and statistical tests",
x = "Plot Fertility Quality",
y = "Average Yield Ratio"
) +
theme_minimal()
###IQR_plot_fert_quality with CA_typology
library(ggplot2)
library(dplyr)
library(scales)
# 1. Clean and prepare data
df_bar <- drivers_data %>%
filter(
!is.na(CA_typology),
!is.na(IQR_plot_fert_quality),
!CA_typology %in% c("no data", "---", "", "NA"),
!IQR_plot_fert_quality %in% c("no data", "---", "", "NA")
) %>%
mutate(
IQR_plot_fert_quality = factor(IQR_plot_fert_quality,
levels = c("high", "medium", "low_quality")) # Desired order
)
# 2. Calculate N per category
n_labels <- df_bar %>%
group_by(IQR_plot_fert_quality) %>%
summarise(n = n(), .groups = "drop")
# 3. Plot
ggplot(df_bar, aes(x = IQR_plot_fert_quality, fill = CA_typology)) +
geom_bar(position = "fill", color = "black") +
geom_text(
data = n_labels,
aes(x = IQR_plot_fert_quality, y = 1.02, label = paste0("N = ", n)),
inherit.aes = FALSE,
size = 4
) +
scale_y_continuous(labels = percent_format(), expand = expansion(mult = c(0, 0.1))) +
labs(
title = "CA_typology distribution by Fertilizer Quality",
subtitle = "Proportions with total sample size per group",
x = "IQR_plot_fert_quality",
y = "Proportion of CA Typology",
fill = "CA Typology"
) +
theme_minimal()
# 4. Chi-squared test
tab <- table(df_bar$IQR_plot_fert_quality, df_bar$CA_typology)
chi_test <- chisq.test(tab)
chi_test##
## Pearson's Chi-squared test
##
## data: tab
## X-squared = 3.8727, df = 6, p-value = 0.6939how does FErtility perceibed by farmers relate to researchers’ indicators of fertility
library(dplyr)
library(ggplot2)
fert_data <- drivers_data %>%
filter(!is.na(IQR_plot_fert_quality)) %>%
mutate(
IQR_plot_fert_quality = droplevels(factor(IQR_plot_fert_quality))
)
## Paste the env_index
library(dplyr)
env_block_avg <- CONV_CA_Ratio %>%
filter(!is.na(IQR_block),
!is.na(env_index)) %>%
group_by(IQR_block) %>%
summarise(
env_index_avg = mean(env_index, na.rm = TRUE),
.groups = "drop"
)
fert_data <- drivers_data %>%
filter(!is.na(IQR_plot_fert_quality)) %>%
mutate(
IQR_plot_fert_quality = droplevels(factor(IQR_plot_fert_quality))
) %>%
left_join(env_block_avg, by = "IQR_block")
cat("Rows in fert_data:", nrow(fert_data), "\n")## Rows in fert_data: 471summary(fert_data$env_index_avg)## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 0.3312 0.8015 0.9828 0.9967 1.1471 2.4047# Check variables relation
# Environmetal index
library(dplyr)
library(ggplot2)
# ---------------------------------------
# 1️⃣ Ensure env_index_avg exists
# ---------------------------------------
fert_data <- fert_data %>%
filter(!is.na(IQR_plot_fert_quality),
!is.na(env_index_avg))
# ---------------------------------------
# 2️⃣ Summary statistics
# ---------------------------------------
summary_env <- fert_data %>%
group_by(IQR_plot_fert_quality) %>%
summarise(
mean_val = mean(env_index_avg, na.rm = TRUE),
N = n(),
y_pos = max(env_index_avg, na.rm = TRUE) * 1.05,
label = paste0("Mean = ", round(mean_val, 2),
"\nN = ", N),
.groups = "drop"
)
# ---------------------------------------
# 3️⃣ Plot
# ---------------------------------------
ggplot(fert_data,
aes(x = IQR_plot_fert_quality,
y = env_index_avg,
fill = IQR_plot_fert_quality)) +
geom_boxplot(alpha = 0.6) +
geom_jitter(width = 0.15,
alpha = 0.3,
color = "black") +
geom_text(data = summary_env,
aes(x = IQR_plot_fert_quality,
y = y_pos,
label = label),
inherit.aes = FALSE,
size = 4) +
labs(
title = "Environmental Index vs Farmer Perceived Soil Fertility",
x = "Farmer Perceived Soil Fertility",
y = "Average Environmental Index (across seasons)"
) +
theme_minimal(base_size = 14) +
theme(
legend.position = "none",
plot.title = element_text(face = "bold", hjust = 0.5)
)
# SOC
library(dplyr)
library(ggplot2)
fert_data <- drivers_data %>%
filter(!is.na(IQR_plot_fert_quality),
!is.na(ICR_Org_C_avg))
summary_OC <- fert_data %>%
group_by(IQR_plot_fert_quality) %>%
summarise(
mean_val = mean(ICR_Org_C_avg, na.rm = TRUE),
N = n(),
y_pos = max(ICR_Org_C_avg, na.rm = TRUE) * 1.05,
label = paste0("Mean=", round(mean_val, 2),
"\nN=", N),
.groups = "drop"
)
# Plot
ggplot(fert_data,
aes(x = IQR_plot_fert_quality,
y = ICR_Org_C_avg,
fill = IQR_plot_fert_quality)) +
geom_boxplot(alpha = 0.6) +
geom_jitter(width = 0.15, alpha = 0.3) +
geom_text(data = summary_OC,
aes(x = IQR_plot_fert_quality,
y = y_pos,
label = label),
inherit.aes = FALSE,
size = 4) +
labs(
title = "Organic Carbon vs Farmer Perception",
x = "Farmer Perceived Soil Fertility",
y = "Organic Carbon"
) +
theme_minimal() +
theme(legend.position = "none")
##P
# Summary
summary_P <- fert_data %>%
group_by(IQR_plot_fert_quality) %>%
summarise(
mean_val = mean(ICR_Avail_P_avg, na.rm = TRUE),
N = n(),
y_pos = max(ICR_Avail_P_avg, na.rm = TRUE) * 1.05,
label = paste0("Mean=", round(mean_val, 1),
"\nN=", N),
.groups = "drop"
)
ggplot(fert_data,
aes(x = IQR_plot_fert_quality,
y = ICR_Avail_P_avg,
fill = IQR_plot_fert_quality)) +
geom_boxplot(alpha = 0.6) +
geom_jitter(width = 0.15, alpha = 0.3) +
geom_text(data = summary_P,
aes(x = IQR_plot_fert_quality,
y = y_pos,
label = label),
inherit.aes = FALSE,
size = 4) +
theme_minimal() +
theme(legend.position = "none")
##K
summary_K <- fert_data %>%
group_by(IQR_plot_fert_quality) %>%
summarise(
mean_val = mean(ICR_K_perc_23A, na.rm = TRUE),
N = n(),
y_pos = max(ICR_K_perc_23A, na.rm = TRUE) * 1.05,
label = paste0("Mean = ", round(mean_val, 3),
"\nN = ", N),
.groups = "drop"
)
# ---------------------------------------
# 3️⃣ Plot
# ---------------------------------------
ggplot(fert_data,
aes(x = IQR_plot_fert_quality,
y = ICR_K_perc_23A,
fill = IQR_plot_fert_quality)) +
geom_boxplot(alpha = 0.6) +
geom_jitter(width = 0.15,
alpha = 0.3,
color = "black") +
geom_text(data = summary_K,
aes(x = IQR_plot_fert_quality,
y = y_pos,
label = label),
inherit.aes = FALSE,
size = 4) +
labs(
title = "Potassium (ICR_K_perc_23A) vs Farmer Perceived Soil Fertility",
x = "Farmer Perceived Soil Fertility",
y = "Potassium (%)"
) +
theme_minimal(base_size = 14) +
theme(
legend.position = "none",
plot.title = element_text(face = "bold", hjust = 0.5)
)
# Soil texture
ggplot(fert_data,
aes(x = IQR_plot_fert_quality,
fill = IQR_soil_texture)) +
geom_bar(position = "fill") +
labs(
title = "Soil Texture Distribution by Farmer Perception",
x = "Farmer Perception",
y = "Proportion",
fill = "Soil Texture"
) +
theme_minimal()
chisq.test(table(fert_data$IQR_plot_fert_quality,
fert_data$IQR_soil_texture))##
## Pearson's Chi-squared test
##
## data: table(fert_data$IQR_plot_fert_quality, fert_data$IQR_soil_texture)
## X-squared = 5.8932, df = 6, p-value = 0.4353# Soil color
ggplot(fert_data,
aes(x = IQR_plot_fert_quality,
fill = IQR_soil_color)) +
geom_bar(position = "fill") +
labs(
title = "Soil Color by Farmer Perception",
y = "Proportion"
) +
theme_minimal()
# Slope position
ggplot(fert_data,
aes(x = IQR_plot_fert_quality,
fill = IQF_position_slope)) +
geom_bar(position = "fill") +
labs(
title = "Slope Position by Farmer Perception",
y = "Proportion"
) +
theme_minimal()
## CA Performance as affected by IQR_soil_texture NOtes: Conclusions:
###IQR_soil_texture with Y_ratio
library(ggplot2)
library(dplyr)
library(multcompView)
# 1. Filter and clean the data
df_filt <- drivers_data %>%
filter(
!is.na(avg_Y_ratio),
!is.na(IQR_soil_texture),
!IQR_soil_texture %in% c("no data", "---", "", "NA")
) %>%
mutate(IQR_soil_texture = droplevels(factor(IQR_soil_texture)))
# Optional: Reorder texture levels manually (if desired)
# Example order: "sand", "loam", "clay", etc.
df_filt <- df_filt %>%
mutate(IQR_soil_texture = factor(IQR_soil_texture, levels = c("Fine", "Mix", "Coarse")))
# 2. Run ANOVA
anova_model <- aov(avg_Y_ratio ~ IQR_soil_texture, data = df_filt)
anova_p <- summary(anova_model)[[1]][["Pr(>F)"]][1]
# 3. Run Kruskal-Wallis test
kruskal_p <- kruskal.test(avg_Y_ratio ~ IQR_soil_texture, data = df_filt)$p.value
# 4. Tukey HSD and compact letter display
tukey <- TukeyHSD(anova_model)
letters <- multcompLetters(tukey$`IQR_soil_texture`[, "p adj"])$Letters
letters_df <- data.frame(IQR_soil_texture = names(letters), Letters = letters)
# 5. Compute group means
means_df <- df_filt %>%
group_by(IQR_soil_texture) %>%
summarise(mean_val = mean(avg_Y_ratio, na.rm = TRUE), .groups = "drop") %>%
mutate(mean_label = paste0("mean = ", round(mean_val, 2)))
# 6. Combine means + letters
labels_df <- left_join(means_df, letters_df, by = "IQR_soil_texture")
# 7. Test label
test_label <- paste0(
"ANOVA p = ", signif(anova_p, 3),
"\nKruskal p = ", signif(kruskal_p, 3)
)
# 8. Plot
ggplot(df_filt, aes(x = IQR_soil_texture, y = avg_Y_ratio)) +
geom_boxplot(outlier.shape = NA, fill = "lightgray") +
geom_jitter(width = 0.2, alpha = 0.6, color = "steelblue") +
# Letters above box
geom_text(
data = labels_df,
aes(x = IQR_soil_texture, y = max(df_filt$avg_Y_ratio, na.rm = TRUE) * 1.02, label = Letters),
inherit.aes = FALSE,
fontface = "bold", size = 5, color = "darkred"
) +
# Means below box
geom_text(
data = labels_df,
aes(x = IQR_soil_texture, y = min(df_filt$avg_Y_ratio, na.rm = TRUE) * 0.95, label = mean_label),
inherit.aes = FALSE,
size = 4
) +
# Statistical tests label
annotate("text",
x = Inf, y = Inf,
hjust = 1.1, vjust = 1.5,
label = test_label,
size = 5) +
labs(
title = "avg_Y_ratio by IQR_soil_texture",
subtitle = "Tukey HSD letters, group means, ANOVA and Kruskal–Wallis p-values",
x = "Soil Texture",
y = "Average Yield Ratio"
) +
theme_minimal()
###IQR_soil_texture with CA_typology
library(ggplot2)
library(dplyr)
library(scales)
# 1. Clean and prepare data
df_bar <- drivers_data %>%
filter(
!is.na(CA_typology),
!is.na(IQR_soil_texture),
!CA_typology %in% c("no data", "---", "", "NA"),
!IQR_soil_texture %in% c("no data", "---", "", "NA")
) %>%
mutate(IQR_soil_texture = droplevels(as.factor(IQR_soil_texture))) %>%
mutate(IQR_soil_texture = factor(IQR_soil_texture ,
levels = c("Fine", "Mix", "Coarse")))
# Desired order
# 2. Calculate N per level of explanatory variable
n_labels <- df_bar %>%
group_by(IQR_soil_texture) %>%
summarise(n = n(), .groups = "drop")
# 3. Bar plot
ggplot(df_bar, aes(x = IQR_soil_texture, fill = CA_typology)) +
geom_bar(position = "fill", color = "black") +
geom_text(
data = n_labels,
aes(x = IQR_soil_texture, y = 1.02, label = paste0("N = ", n)),
inherit.aes = FALSE,
size = 4
) +
scale_y_continuous(labels = percent_format(), expand = expansion(mult = c(0, 0.1))) +
labs(
title = "CA_typology distribution by Soil Texture",
subtitle = "Proportions with total sample size per texture class",
x = "IQR_soil_texture",
y = "Proportion of CA Typology",
fill = "CA Typology"
) +
theme_minimal()
# 4. Statistical tests
tab <- table(df_bar$IQR_soil_texture, df_bar$CA_typology)
# Chi-squared test
chi_test <- chisq.test(tab)
chi_test##
## Pearson's Chi-squared test
##
## data: tab
## X-squared = 6.0014, df = 6, p-value = 0.423CA Performance as affected by IQR_soil_color
NOtes: Conclusions: ###IQR_soil_color with Y_ratio
library(ggplot2)
library(dplyr)
library(multcompView)
# 1. Filter and clean the data
df_filt <- drivers_data %>%
filter(
!is.na(avg_Y_ratio),
!is.na(IQR_soil_color),
!IQR_soil_color %in% c("no data", "---", "", "#N/A")
) %>%
mutate(IQR_soil_color = droplevels(factor(IQR_soil_color)))
# 2. Run ANOVA
anova_model <- aov(avg_Y_ratio ~ IQR_soil_color, data = df_filt)
anova_p <- summary(anova_model)[[1]][["Pr(>F)"]][1]
# 3. Kruskal–Wallis test
kruskal_p <- kruskal.test(avg_Y_ratio ~ IQR_soil_color, data = df_filt)$p.value
# 4. Tukey HSD and letters
tukey <- TukeyHSD(anova_model)
letters <- multcompLetters(tukey$`IQR_soil_color`[, "p adj"])$Letters
letters_df <- data.frame(IQR_soil_color = names(letters), Letters = letters)
# 5. Group means
means_df <- df_filt %>%
group_by(IQR_soil_color) %>%
summarise(mean_val = mean(avg_Y_ratio, na.rm = TRUE), .groups = "drop") %>%
mutate(mean_label = paste0("mean = ", round(mean_val, 2)))
# 6. Combine letters and means
labels_df <- left_join(means_df, letters_df, by = "IQR_soil_color")
# 7. Test label
test_label <- paste0(
"ANOVA p = ", signif(anova_p, 3),
"\nKruskal p = ", signif(kruskal_p, 3)
)
# 8. Plot
ggplot(df_filt, aes(x = IQR_soil_color, y = avg_Y_ratio)) +
geom_boxplot(outlier.shape = NA, fill = "lightgray") +
geom_jitter(width = 0.2, alpha = 0.6, color = "steelblue") +
# Letters above boxes
geom_text(
data = labels_df,
aes(x = IQR_soil_color, y = max(df_filt$avg_Y_ratio, na.rm = TRUE) * 1.02, label = Letters),
inherit.aes = FALSE,
fontface = "bold", size = 5, color = "darkred"
) +
# Means below boxes
geom_text(
data = labels_df,
aes(x = IQR_soil_color, y = min(df_filt$avg_Y_ratio, na.rm = TRUE) * 0.95, label = mean_label),
inherit.aes = FALSE,
size = 4
) +
annotate("text",
x = Inf, y = Inf,
hjust = 1.1, vjust = 1.5,
label = test_label,
size = 5) +
labs(
title = "avg_Y_ratio by IQR_soil_color",
subtitle = "Tukey HSD letters, group means, ANOVA and Kruskal–Wallis p-values",
x = "Soil Color",
y = "Average Yield Ratio"
) +
theme_minimal()
### Soil color facet by environment
library(ggplot2)
library(dplyr)
library(multcompView)
# -------------------------------------------------
# 1️⃣ Filter and clean
# -------------------------------------------------
df_filt <- drivers_data %>%
filter(
!is.na(avg_Y_ratio),
!is.na(IQR_soil_color),
!is.na(IQF_environment),
!IQR_soil_color %in% c("no data", "---", "", "#N/A")
) %>%
mutate(IQR_soil_color = droplevels(factor(IQR_soil_color)))
# -------------------------------------------------
# 2️⃣ Run ANOVA + Kruskal + Tukey per environment
# -------------------------------------------------
results_env <- df_filt %>%
group_by(IQF_environment) %>%
group_modify(~{
# Count number of soil color groups
n_groups <- length(unique(.x$IQR_soil_color))
# Compute group means (always)
means_df <- .x %>%
group_by(IQR_soil_color) %>%
summarise(
mean_val = mean(avg_Y_ratio, na.rm = TRUE),
.groups = "drop"
) %>%
mutate(mean_label = paste0("mean = ", round(mean_val, 2)))
# If fewer than 2 groups → skip tests
if (n_groups < 2) {
means_df$Letters <- ""
means_df$anova_p <- NA
means_df$kruskal_p <- NA
return(means_df)
}
# ---------------- ANOVA ----------------
anova_model <- aov(avg_Y_ratio ~ IQR_soil_color, data = .x)
anova_p <- summary(anova_model)[[1]][["Pr(>F)"]][1]
# ---------------- Kruskal ----------------
kruskal_p <- kruskal.test(avg_Y_ratio ~ IQR_soil_color, data = .x)$p.value
# ---------------- Tukey (robust) ----------------
tukey <- TukeyHSD(anova_model)
tukey_mat <- tukey$IQR_soil_color
if (is.null(tukey_mat) || nrow(tukey_mat) == 0) {
means_df$Letters <- ""
} else {
tukey_vec <- tukey_mat[, "p adj"]
# Ensure names exist
if (is.null(names(tukey_vec))) {
names(tukey_vec) <- rownames(tukey_mat)
}
letters <- multcompLetters(tukey_vec)$Letters
letters_df <- data.frame(
IQR_soil_color = names(letters),
Letters = letters,
stringsAsFactors = FALSE
)
means_df <- left_join(means_df, letters_df,
by = "IQR_soil_color")
}
means_df$anova_p <- anova_p
means_df$kruskal_p <- kruskal_p
means_df
}) %>%
ungroup()
# -------------------------------------------------
# 3️⃣ Create p-value labels per environment
# -------------------------------------------------
test_labels <- results_env %>%
distinct(IQF_environment, anova_p, kruskal_p) %>%
mutate(
test_label = paste0(
"ANOVA p = ", signif(anova_p, 3),
"\nKruskal p = ", signif(kruskal_p, 3)
)
)
# -------------------------------------------------
# 4️⃣ Plot
# -------------------------------------------------
ggplot(df_filt, aes(x = IQR_soil_color, y = avg_Y_ratio)) +
geom_boxplot(outlier.shape = NA, fill = "lightgray") +
geom_jitter(width = 0.2, alpha = 0.6, color = "steelblue") +
# Tukey Letters
geom_text(
data = results_env,
aes(
x = IQR_soil_color,
y = max(df_filt$avg_Y_ratio, na.rm = TRUE) * 1.05,
label = Letters
),
inherit.aes = FALSE,
fontface = "bold",
size = 5,
color = "darkred"
) +
# Means
geom_text(
data = results_env,
aes(
x = IQR_soil_color,
y = min(df_filt$avg_Y_ratio, na.rm = TRUE) * 0.95,
label = mean_label
),
inherit.aes = FALSE,
size = 4
) +
# P-values
geom_text(
data = test_labels,
aes(x = Inf, y = Inf, label = test_label),
hjust = 1.1,
vjust = 1.5,
inherit.aes = FALSE,
size = 4
) +
facet_wrap(~ IQF_environment) +
labs(
title = "avg_Y_ratio by IQR_soil_color",
subtitle = "Facet by Environment | Tukey letters + ANOVA + Kruskal–Wallis",
x = "Soil Color",
y = "Average Yield Ratio"
) +
theme_minimal()
###IQR_soil_color with CA_typology
library(ggplot2)
library(dplyr)
library(scales)
# 1. Clean data
df_bar <- drivers_data %>%
filter(
!is.na(CA_typology),
!is.na(IQR_soil_color),
!CA_typology %in% c("no data", "---", "", "NA"),
!IQR_soil_color %in% c("no data", "---", "", "#N/A")
)
# 2. Calculate N per category
n_labels <- df_bar %>%
group_by(IQR_soil_color) %>%
summarise(n = n(), .groups = "drop")
# 3. Plot
ggplot(df_bar, aes(x = IQR_soil_color, fill = CA_typology)) +
geom_bar(position = "fill", color = "black") +
geom_text(
data = n_labels,
aes(x = IQR_soil_color, y = 1.02, label = paste0("N = ", n)),
inherit.aes = FALSE,
size = 4
) +
scale_y_continuous(labels = percent_format(), expand = expansion(mult = c(0, 0.1))) +
labs(
title = "CA_typology distribution by Soil Color",
subtitle = "Proportional bars with total sample size",
x = "Soil Color",
y = "Proportion of CA Typology",
fill = "CA Typology"
) +
theme_minimal()
# 4. Chi-squared and Fisher’s test
tab <- table(df_bar$IQR_soil_color, df_bar$CA_typology)
# Chi-squared test
chi_test <- chisq.test(tab)
print(chi_test)##
## Pearson's Chi-squared test
##
## data: tab
## X-squared = 6.7453, df = 9, p-value = 0.6636# Fisher’s Exact Test (with simulation to avoid memory issues)
fisher_test <- fisher.test(tab, simulate.p.value = TRUE, B = 1e5)
print(fisher_test)##
## Fisher's Exact Test for Count Data with simulated p-value (based on
## 1e+05 replicates)
##
## data: tab
## p-value = 0.6509
## alternative hypothesis: two.sidedCA Performance as affected by IQF_position_slope
NOtes: Conclusions: ###IQF_position_slope with Y_ratio
library(ggplot2)
library(dplyr)
library(multcompView)
# 1. Filter and clean data
df_filt <- drivers_data %>%
filter(
!is.na(avg_Y_ratio),
!is.na(IQF_position_slope),
!IQF_position_slope %in% c("no data", "---", "", "NA")
) %>%
mutate(IQF_position_slope = droplevels(factor(IQF_position_slope)))
# 2. ANOVA
anova_model <- aov(avg_Y_ratio ~ IQF_position_slope, data = df_filt)
anova_p <- summary(anova_model)[[1]][["Pr(>F)"]][1]
# 3. Kruskal-Wallis test
kruskal_p <- kruskal.test(avg_Y_ratio ~ IQF_position_slope, data = df_filt)$p.value
# 4. Tukey HSD and group letters
tukey <- TukeyHSD(anova_model)
letters <- multcompLetters(tukey$`IQF_position_slope`[, "p adj"])$Letters
letters_df <- data.frame(IQF_position_slope = names(letters), Letters = letters)
# 5. Group means
means_df <- df_filt %>%
group_by(IQF_position_slope) %>%
summarise(mean_val = mean(avg_Y_ratio, na.rm = TRUE), .groups = "drop") %>%
mutate(mean_label = paste0("mean = ", round(mean_val, 2)))
# 6. Combine means + letters
labels_df <- left_join(means_df, letters_df, by = "IQF_position_slope")
# 7. Statistical test label
test_label <- paste0(
"ANOVA p = ", signif(anova_p, 3),
"\nKruskal p = ", signif(kruskal_p, 3)
)
# 8. Plot
ggplot(df_filt, aes(x = IQF_position_slope, y = avg_Y_ratio)) +
geom_boxplot(outlier.shape = NA, fill = "lightgray") +
geom_jitter(width = 0.2, alpha = 0.6, color = "steelblue") +
# Tukey HSD letters (above)
geom_text(
data = labels_df,
aes(x = IQF_position_slope, y = max(df_filt$avg_Y_ratio, na.rm = TRUE) * 1.02, label = Letters),
inherit.aes = FALSE,
fontface = "bold", size = 5, color = "darkred"
) +
# Means (below)
geom_text(
data = labels_df,
aes(x = IQF_position_slope, y = min(df_filt$avg_Y_ratio, na.rm = TRUE) * 0.95, label = mean_label),
inherit.aes = FALSE,
size = 4
) +
annotate("text",
x = Inf, y = Inf,
hjust = 1.1, vjust = 1.5,
label = test_label,
size = 5) +
labs(
title = "avg_Y_ratio by IQF_position_slope",
subtitle = "Tukey HSD letters, group means, ANOVA and Kruskal–Wallis p-values",
x = "Slope Position",
y = "Average Yield Ratio"
) +
theme_minimal()
###Facet by environment
library(ggplot2)
library(dplyr)
library(multcompView)
# -------------------------------------------------
# 1️⃣ Filter and clean
# -------------------------------------------------
df_filt <- drivers_data %>%
filter(
!is.na(avg_Y_ratio),
!is.na(IQF_position_slope),
!is.na(IQF_environment),
!IQF_position_slope %in% c("no data", "---", "", "NA")
) %>%
mutate(IQF_position_slope = droplevels(factor(IQF_position_slope)))
# -------------------------------------------------
# 2️⃣ Run ANOVA + Kruskal + Tukey per environment
# -------------------------------------------------
results_env <- df_filt %>%
group_by(IQF_environment) %>%
group_modify(~{
n_groups <- length(unique(.x$IQF_position_slope))
# Means
means_df <- .x %>%
group_by(IQF_position_slope) %>%
summarise(
mean_val = mean(avg_Y_ratio, na.rm = TRUE),
.groups = "drop"
) %>%
mutate(mean_label = paste0("mean = ", round(mean_val, 2)))
if (n_groups < 2) {
means_df$Letters <- ""
means_df$anova_p <- NA
means_df$kruskal_p <- NA
return(means_df)
}
# ANOVA
anova_model <- aov(avg_Y_ratio ~ IQF_position_slope, data = .x)
anova_p <- summary(anova_model)[[1]][["Pr(>F)"]][1]
# Kruskal
kruskal_p <- kruskal.test(avg_Y_ratio ~ IQF_position_slope, data = .x)$p.value
# Tukey (robust)
tukey <- TukeyHSD(anova_model)
tukey_mat <- tukey$IQF_position_slope
if (is.null(tukey_mat) || nrow(tukey_mat) == 0) {
means_df$Letters <- ""
} else {
tukey_vec <- tukey_mat[, "p adj"]
if (is.null(names(tukey_vec))) {
names(tukey_vec) <- rownames(tukey_mat)
}
letters <- multcompLetters(tukey_vec)$Letters
letters_df <- data.frame(
IQF_position_slope = names(letters),
Letters = letters,
stringsAsFactors = FALSE
)
means_df <- left_join(means_df, letters_df,
by = "IQF_position_slope")
}
means_df$anova_p <- anova_p
means_df$kruskal_p <- kruskal_p
means_df
}) %>%
ungroup()
# -------------------------------------------------
# 3️⃣ Create p-value labels per environment
# -------------------------------------------------
test_labels <- results_env %>%
distinct(IQF_environment, anova_p, kruskal_p) %>%
mutate(
test_label = paste0(
"ANOVA p = ", signif(anova_p, 3),
"\nKruskal p = ", signif(kruskal_p, 3)
)
)
# -------------------------------------------------
# 4️⃣ Plot
# -------------------------------------------------
ggplot(df_filt, aes(x = IQF_position_slope, y = avg_Y_ratio)) +
geom_boxplot(outlier.shape = NA, fill = "lightgray") +
geom_jitter(width = 0.2, alpha = 0.6, color = "steelblue") +
# Tukey Letters
geom_text(
data = results_env,
aes(
x = IQF_position_slope,
y = max(df_filt$avg_Y_ratio, na.rm = TRUE) * 1.05,
label = Letters
),
inherit.aes = FALSE,
fontface = "bold",
size = 5,
color = "darkred"
) +
# Means
geom_text(
data = results_env,
aes(
x = IQF_position_slope,
y = min(df_filt$avg_Y_ratio, na.rm = TRUE) * 0.95,
label = mean_label
),
inherit.aes = FALSE,
size = 4
) +
# P-values
geom_text(
data = test_labels,
aes(x = Inf, y = Inf, label = test_label),
hjust = 1.1,
vjust = 1.5,
inherit.aes = FALSE,
size = 4
) +
facet_wrap(~ IQF_environment) +
labs(
title = "avg_Y_ratio by IQF_position_slope",
subtitle = "Facet by Environment | Tukey letters + ANOVA + Kruskal–Wallis",
x = "Slope Position",
y = "Average Yield Ratio"
) +
theme_minimal()
###IQF_position_slope with CA_typology
library(ggplot2)
library(dplyr)
library(scales)
# 1. Clean and prepare data
df_bar <- drivers_data %>%
filter(
!is.na(CA_typology),
!is.na(IQF_position_slope),
!CA_typology %in% c("no data", "---", "", "NA"),
!IQF_position_slope %in% c("no data", "---", "", "NA")
)
# 2. Calculate total N per slope category
n_labels <- df_bar %>%
group_by(IQF_position_slope) %>%
summarise(n = n(), .groups = "drop")
# 3. Plot proportion bar chart
ggplot(df_bar, aes(x = IQF_position_slope, fill = CA_typology)) +
geom_bar(position = "fill", color = "black") +
geom_text(
data = n_labels,
aes(x = IQF_position_slope, y = 1.02, label = paste0("N = ", n)),
inherit.aes = FALSE,
size = 4
) +
scale_y_continuous(labels = percent_format(), expand = expansion(mult = c(0, 0.1))) +
labs(
title = "CA_typology distribution by Position on Slope",
subtitle = "Proportions and total N per slope category",
x = "Slope Position",
y = "Proportion of CA Typology",
fill = "CA Typology"
) +
theme_minimal()
# 4. Run statistical tests
tab <- table(df_bar$IQF_position_slope, df_bar$CA_typology)
# Chi-squared test
chi_test <- chisq.test(tab)
print(chi_test)##
## Pearson's Chi-squared test
##
## data: tab
## X-squared = 19.612, df = 12, p-value = 0.07478# Fisher’s Exact Test with simulation fallback
fisher_test <- fisher.test(tab, simulate.p.value = TRUE, B = 1e5)
print(fisher_test)##
## Fisher's Exact Test for Count Data with simulated p-value (based on
## 1e+05 replicates)
##
## data: tab
## p-value = 0.1505
## alternative hypothesis: two.sidedCA Performance as affected by Weed_species_A_combined
NOtes: Conclusions: ###Weed_species_A_combined with Y_ratio
library(ggplot2)
library(dplyr)
library(multcompView)
# 1. Filter and clean data
df_filt <- drivers_data %>%
filter(
!is.na(avg_Y_ratio),
!is.na(Weed_species_A_combined),
!Weed_species_A_combined %in% c("no data", "-", "—", "NA")
) %>%
mutate(Weed_species_A_combined = droplevels(factor(Weed_species_A_combined)))
# 2. ANOVA
anova_model <- aov(avg_Y_ratio ~ Weed_species_A_combined, data = df_filt)
anova_p <- summary(anova_model)[[1]][["Pr(>F)"]][1]
# 3. Kruskal–Wallis
kruskal_p <- kruskal.test(avg_Y_ratio ~ Weed_species_A_combined, data = df_filt)$p.value
# 4. Tukey HSD and group letters
tukey <- TukeyHSD(anova_model)
letters <- multcompLetters(tukey$`Weed_species_A_combined`[, "p adj"])$Letters
letters_df <- data.frame(Weed_species_A_combined = names(letters), Letters = letters)
# 5. Group means
means_df <- df_filt %>%
group_by(Weed_species_A_combined) %>%
summarise(mean_val = mean(avg_Y_ratio, na.rm = TRUE), .groups = "drop") %>%
mutate(mean_label = paste0("mean = ", round(mean_val, 2)))
# 6. Combine means + letters
labels_df <- left_join(means_df, letters_df, by = "Weed_species_A_combined")
# 7. Combine test labels
test_label <- paste0(
"ANOVA p = ", signif(anova_p, 3),
"\nKruskal p = ", signif(kruskal_p, 3)
)
# 8. Plot
ggplot(df_filt, aes(x = Weed_species_A_combined, y = avg_Y_ratio)) +
geom_boxplot(outlier.shape = NA, fill = "lightgray") +
geom_jitter(width = 0.2, alpha = 0.6, color = "steelblue") +
# Letters above
geom_text(
data = labels_df,
aes(x = Weed_species_A_combined, y = max(df_filt$avg_Y_ratio, na.rm = TRUE) * 1.02, label = Letters),
inherit.aes = FALSE,
fontface = "bold", size = 5, color = "darkred"
) +
# Means below
geom_text(
data = labels_df,
aes(x = Weed_species_A_combined, y = min(df_filt$avg_Y_ratio, na.rm = TRUE) * 0.95, label = mean_label),
inherit.aes = FALSE,
size = 4
) +
annotate("text",
x = Inf, y = Inf,
hjust = 1.1, vjust = 1.5,
label = test_label,
size = 5) +
labs(
title = "avg_Y_ratio by Weed_species_A_combined",
subtitle = "Tukey HSD letters, group means, ANOVA and Kruskal–Wallis p-values",
x = "Weed Species (Combined)",
y = "Average Yield Ratio"
) +
theme_minimal()
###Weed_species_A_combined with CA_typology
library(ggplot2)
library(dplyr)
library(scales)
# 1. Filter and clean data
df_bar <- drivers_data %>%
filter(
!is.na(CA_typology),
!is.na(Weed_species_A_combined),
!CA_typology %in% c("no data", "---", "", "NA"),
!Weed_species_A_combined %in% c("no data", "—", "", "NA")
)
# 2. Total N per weed category
n_labels <- df_bar %>%
group_by(Weed_species_A_combined) %>%
summarise(n = n(), .groups = "drop")
# 3. Proportional bar plot
ggplot(df_bar, aes(x = Weed_species_A_combined, fill = CA_typology)) +
geom_bar(position = "fill", color = "black") +
geom_text(
data = n_labels,
aes(x = Weed_species_A_combined, y = 1.02, label = paste0("N = ", n)),
inherit.aes = FALSE,
size = 4
) +
scale_y_continuous(labels = percent_format(), expand = expansion(mult = c(0, 0.1))) +
labs(
title = "CA_typology distribution by Dominant Weed Species",
subtitle = "Stacked bars show % CA_typology, with N on top",
x = "Weed Species Group (A combined)",
y = "Proportion of CA Typology",
fill = "CA Typology"
) +
theme_minimal()
# 4. Statistical tests
tab <- table(df_bar$Weed_species_A_combined, df_bar$CA_typology)
# Chi-squared test
chi_test <- chisq.test(tab)
print(chi_test)##
## Pearson's Chi-squared test
##
## data: tab
## X-squared = 5.9079, df = 6, p-value = 0.4336# Fisher’s Exact Test with simulated p-value (for large tables)
fisher_test <- fisher.test(tab, simulate.p.value = TRUE, B = 1e5)
print(fisher_test)##
## Fisher's Exact Test for Count Data with simulated p-value (based on
## 1e+05 replicates)
##
## data: tab
## p-value = 0.3731
## alternative hypothesis: two.sidedCA Performance as affected by IQF_environment
NOtes: Conclusions: ###IQF_environment with Y_ratio
library(ggplot2)
library(dplyr)
library(multcompView)
# 1. Filter and clean data
df_filt <- drivers_data %>%
filter(
!is.na(avg_Y_ratio),
!is.na(IQF_environment),
!IQF_environment %in% c("no data", "---", "", "NA")
) %>%
mutate(IQF_environment = droplevels(factor(IQF_environment)))
df_filt <- df_filt %>%
mutate(IQF_environment = factor(IQF_environment, levels = c("Highland", "Transition", "Mid_Low_Wet", "Mid_Low_Dry")))
# 2. ANOVA
anova_model <- aov(avg_Y_ratio ~ IQF_environment, data = df_filt)
anova_p <- summary(anova_model)[[1]][["Pr(>F)"]][1]
# 3. Kruskal-Wallis test
kruskal_p <- kruskal.test(avg_Y_ratio ~ IQF_environment, data = df_filt)$p.value
# 4. Tukey HSD + Letters
tukey <- TukeyHSD(anova_model)
letters <- multcompLetters(tukey$`IQF_environment`[, "p adj"])$Letters
letters_df <- data.frame(IQF_environment = names(letters), Letters = letters)
# 5. Group means
means_df <- df_filt %>%
group_by(IQF_environment) %>%
summarise(mean_val = mean(avg_Y_ratio, na.rm = TRUE), .groups = "drop") %>%
mutate(mean_label = paste0("mean = ", round(mean_val, 2)))
# 6. Combine letters + means
labels_df <- left_join(means_df, letters_df, by = "IQF_environment")
# 7. Test summary
test_label <- paste0(
"ANOVA p = ", signif(anova_p, 3),
"\nKruskal p = ", signif(kruskal_p, 3)
)
# 8. Plot
ggplot(df_filt, aes(x = IQF_environment, y = avg_Y_ratio)) +
geom_boxplot(outlier.shape = NA, fill = "lightgray") +
geom_jitter(width = 0.2, alpha = 0.6, color = "steelblue") +
# Letters above
geom_text(
data = labels_df,
aes(x = IQF_environment, y = max(df_filt$avg_Y_ratio, na.rm = TRUE) * 1.02, label = Letters),
inherit.aes = FALSE,
fontface = "bold", size = 5, color = "darkred"
) +
# Means below
geom_text(
data = labels_df,
aes(x = IQF_environment, y = min(df_filt$avg_Y_ratio, na.rm = TRUE) * 0.95, label = mean_label),
inherit.aes = FALSE,
size = 4
) +
annotate("text",
x = Inf, y = Inf,
hjust = 1.1, vjust = 1.5,
label = test_label,
size = 5) +
labs(
title = "avg_Y_ratio by IQF_environment",
subtitle = "Tukey HSD letters, group means, ANOVA and Kruskal–Wallis p-values",
x = "Environment Type",
y = "Average Yield Ratio"
) +
theme_minimal()
###IQF_environment with CA_typology
library(ggplot2)
library(dplyr)
library(scales)
# 1. Clean data
df_bar <- drivers_data %>%
filter(
!is.na(CA_typology),
!is.na(IQF_environment),
!CA_typology %in% c("no data", "---", "", "NA"),
!IQF_environment %in% c("no data", "---", "", "NA")
)
df_bar <- df_bar %>%
mutate(IQF_environment = factor(IQF_environment, levels = c("Highland", "Transition", "Mid_Low_Wet", "Mid_Low_Dry")))
# 2. Total N per environment
n_labels <- df_bar %>%
group_by(IQF_environment) %>%
summarise(n = n(), .groups = "drop")
# 3. Proportional bar plot
ggplot(df_bar, aes(x = IQF_environment, fill = CA_typology)) +
geom_bar(position = "fill", color = "black") +
geom_text(
data = n_labels,
aes(x = IQF_environment, y = 1.02, label = paste0("N = ", n)),
inherit.aes = FALSE,
size = 4
) +
scale_y_continuous(labels = percent_format(), expand = expansion(mult = c(0, 0.1))) +
labs(
title = "CA_typology Distribution by Environment",
subtitle = "Proportional bars with total N per environment",
x = "IQF Environment",
y = "Proportion of CA Typology",
fill = "CA Typology"
) +
theme_minimal()
# 4. Statistical tests
tab <- table(df_bar$IQF_environment, df_bar$CA_typology)
# Chi-squared test
chi_test <- chisq.test(tab)
print(chi_test)##
## Pearson's Chi-squared test
##
## data: tab
## X-squared = 25.835, df = 9, p-value = 0.002174# Fisher’s Exact Test with simulation
fisher_test <- fisher.test(tab, simulate.p.value = TRUE, B = 1e5)
print(fisher_test)##
## Fisher's Exact Test for Count Data with simulated p-value (based on
## 1e+05 replicates)
##
## data: tab
## p-value = 0.00789
## alternative hypothesis: two.sidedWe evaluate some variables that where initially dropped
Total residues, Termite insidence, couch grass as weed and Erosion reduction ### Exploration of Total green residues at harvest variable
library(dplyr)
# Step 1: Calculate total unique blocks per season
block_summary <- CONV_CA_Ratio %>%
group_by(IQR_Season) %>%
summarise(
total_blocks = n_distinct(IQR_block),
blocks_with_residues_gt0 = n_distinct(IQR_block[ICR_total_residues_tn_ha > 0]),
blocks_with_residues_nonNA = n_distinct(IQR_block[!is.na(ICR_total_residues_tn_ha)]),
.groups = "drop"
) %>%
mutate(
perc_gt0 = 100 * blocks_with_residues_gt0 / total_blocks,
perc_nonNA = 100 * blocks_with_residues_nonNA / total_blocks
)
# View the result
print(block_summary)## # A tibble: 6 × 6
## IQR_Season total_blocks blocks_with_residues…¹ blocks_with_residues…² perc_gt0
## <chr> <int> <int> <int> <dbl>
## 1 23A 538 537 536 99.8
## 2 23B 499 234 499 46.9
## 3 24A 449 419 449 93.3
## 4 24B 463 461 463 99.6
## 5 25A 428 427 428 99.8
## 6 25B 441 441 441 100
## # ℹ abbreviated names: ¹blocks_with_residues_gt0, ²blocks_with_residues_nonNA
## # ℹ 1 more variable: perc_nonNA <dbl>library(dplyr)
residue_summary <- CONV_CA_Ratio %>%
group_by(IQR_Season) %>%
summarise(
count = sum(!is.na(ICR_total_residues_tn_ha)),
min = min(ICR_total_residues_tn_ha, na.rm = TRUE),
q25 = quantile(ICR_total_residues_tn_ha, 0.25, na.rm = TRUE),
median = median(ICR_total_residues_tn_ha, na.rm = TRUE),
mean = mean(ICR_total_residues_tn_ha, na.rm = TRUE),
q75 = quantile(ICR_total_residues_tn_ha, 0.75, na.rm = TRUE),
max = max(ICR_total_residues_tn_ha, na.rm = TRUE),
sd = sd(ICR_total_residues_tn_ha, na.rm = TRUE)
) %>%
arrange(IQR_Season)
print(residue_summary)## # A tibble: 6 × 9
## IQR_Season count min q25 median mean q75 max sd
## <chr> <int> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
## 1 23A 536 3 11.7 15.9 16.7 20.0 139. 8.37
## 2 23B 499 0 0 0 0.374 0.670 2.82 0.520
## 3 24A 449 0 8 16.3 17.1 24.5 189. 13.8
## 4 24B 463 0 0.6 1.10 1.87 1.86 124. 5.94
## 5 25A 428 0 12.0 20.2 23.4 32.0 141. 15.8
## 6 25B 441 0.0472 0.675 1.06 1.70 1.73 174. 8.27library(ggplot2)
ggplot(CONV_CA_Ratio, aes(x = IQR_Season, y = ICR_total_residues_tn_ha)) +
geom_boxplot(outlier.color = "red", fill = "skyblue") +
labs(title = "Distribution of Total Residues by Season",
y = "ICR_total_residues_tn_ha", x = "IQR_Season") +
theme_minimal()
Exploration of Total perceived mulch at planting variable
library(dplyr)
block_summary <- CONV_CA_Ratio %>%
group_by(IQR_Season) %>%
summarise(
total_blocks = n_distinct(IQR_block),
blocks_with_mulch_gt0 = n_distinct(IQR_block[ICR_mulch_perc_planti > 0]),
blocks_with_mulch_nonNA = n_distinct(IQR_block[!is.na(ICR_mulch_perc_planti)]),
.groups = "drop"
) %>%
mutate(
perc_gt0 = 100 * blocks_with_mulch_gt0 / total_blocks,
perc_nonNA = 100 * blocks_with_mulch_nonNA / total_blocks
)
print(block_summary)## # A tibble: 6 × 6
## IQR_Season total_blocks blocks_with_mulch_gt0 blocks_with_mulch_nonNA perc_gt0
## <chr> <int> <int> <int> <dbl>
## 1 23A 538 0 538 0
## 2 23B 499 479 499 96.0
## 3 24A 449 374 449 83.3
## 4 24B 463 446 463 96.3
## 5 25A 428 308 426 72.0
## 6 25B 441 430 429 97.5
## # ℹ 1 more variable: perc_nonNA <dbl>library(dplyr)
mulch_summary <- CONV_CA_Ratio %>%
group_by(IQR_Season) %>%
summarise(
count = sum(!is.na(ICR_mulch_perc_planti)),
min = min(ICR_mulch_perc_planti, na.rm = TRUE),
q25 = quantile(ICR_mulch_perc_planti, 0.25, na.rm = TRUE),
median = median(ICR_mulch_perc_planti, na.rm = TRUE),
mean = mean(ICR_mulch_perc_planti, na.rm = TRUE),
q75 = quantile(ICR_mulch_perc_planti, 0.75, na.rm = TRUE),
max = max(ICR_mulch_perc_planti, na.rm = TRUE),
sd = sd(ICR_mulch_perc_planti, na.rm = TRUE)
) %>%
arrange(IQR_Season)
print(mulch_summary)## # A tibble: 6 × 9
## IQR_Season count min q25 median mean q75 max sd
## <chr> <int> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
## 1 23A 538 0 0 0 0 0 0 0
## 2 23B 499 0 46.5 65 56.2 75 162 28.2
## 3 24A 449 0 15 30 28.4 40 100 22.4
## 4 24B 463 0 30 50 49.6 70 100 26.9
## 5 25A 426 0 0 20 22.5 35 80 20.3
## 6 25B 429 5 40 60 57.3 75 100 22.0library(ggplot2)
ggplot(CONV_CA_Ratio, aes(x = IQR_Season, y = ICR_mulch_perc_planti)) +
geom_boxplot(outlier.color = "red", fill = "skyblue") +
labs(title = "Distribution of Mulch per Plant by Season",
y = "ICR_mulch_per_planti",
x = "IQR_Season") +
theme_minimal()
model <- lm(Y_ratio ~ ICR_mulch_perc_planti,
data = CONV_CA_Ratio,
na.action = na.omit)
summary(model)##
## Call:
## lm(formula = Y_ratio ~ ICR_mulch_perc_planti, data = CONV_CA_Ratio,
## na.action = na.omit)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.7573 -0.2142 -0.0365 0.1512 3.3184
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 0.8018102 0.0101032 79.362 <2e-16 ***
## ICR_mulch_perc_planti 0.0003911 0.0002177 1.796 0.0725 .
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.3523 on 2802 degrees of freedom
## (14 observations deleted due to missingness)
## Multiple R-squared: 0.00115, Adjusted R-squared: 0.000794
## F-statistic: 3.227 on 1 and 2802 DF, p-value: 0.07252library(ggplot2)
library(ggplot2)
library(dplyr)
library(ggplot2)
data_filt <- CONV_CA_Ratio %>%
filter(year != 23)
ggplot(data_filt,
aes(x = ICR_mulch_perc_planti,
y = Y_ratio,
color = crop)) +
geom_point(alpha = 0.6) +
geom_smooth(method = "lm", se = TRUE) +
labs(
title = "Yield Ratio vs Mulch Percentage per Plant by Crop (Year 23 removed)",
x = "ICR_mulch_perc_planti",
y = "Y_ratio",
color = "Crop"
) +
theme_minimal()
### Is there an optimum much cover? Enough to reach close to plateau NO,
it gous up utill 100% coverage
library(dplyr)
library(ggplot2)
library(segmented)
library(purrr)
# 1️⃣ Clean data (remove year 23 if desired)
data_filt <- CONV_CA_Ratio %>%
filter(year != 23,
!is.na(Y_ratio),
!is.na(ICR_mulch_perc_planti),
!is.na(crop))
# 2️⃣ Fit segmented model per crop
seg_models <- data_filt %>%
group_by(crop) %>%
group_modify(~{
# Linear model first
lm_mod <- lm(Y_ratio ~ ICR_mulch_perc_planti, data = .x)
# Segmented model (starting guess at 30%)
seg_mod <- tryCatch(
segmented(lm_mod,
seg.Z = ~ ICR_mulch_perc_planti,
psi = list(ICR_mulch_perc_planti = 30)),
error = function(e) NULL
)
if (is.null(seg_mod)) {
.x$pred <- predict(lm_mod)
.x$breakpoint <- NA
} else {
.x$pred <- predict(seg_mod)
.x$breakpoint <- seg_mod$psi[2]
}
.x
}) %>%
ungroup()
# 3️⃣ Extract breakpoint per crop
breakpoints <- seg_models %>%
group_by(crop) %>%
summarise(breakpoint = first(breakpoint), .groups = "drop")
print(breakpoints)## # A tibble: 2 × 2
## crop breakpoint
## <chr> <dbl>
## 1 Beans 15.0
## 2 Maize 30.0# 4️⃣ Plot with facet by crop
ggplot(seg_models,
aes(x = ICR_mulch_perc_planti,
y = Y_ratio)) +
geom_point(alpha = 0.4, color = "steelblue") +
geom_line(aes(y = pred),
color = "red",
size = 1.2) +
facet_wrap(~ crop, scales = "free") +
labs(
title = "Yield Ratio vs Mulch Percentage",
subtitle = "Segmented (Plateau) Model per Crop",
x = "Mulch (%)",
y = "Y_ratio"
) +
theme_minimal(base_size = 14)
Mulch ICR_mulch_thicknesscm_planting
library(dplyr)
block_summary <- CONV_CA_Ratio %>%
group_by(IQR_Season) %>%
summarise(
total_blocks = n_distinct(IQR_block),
blocks_with_mulch_gt0 = n_distinct(IQR_block[ICR_mulch_thicknesscm_planting > 0]),
blocks_with_mulch_nonNA = n_distinct(IQR_block[!is.na(ICR_mulch_thicknesscm_planting)]),
.groups = "drop"
) %>%
mutate(
perc_gt0 = 100 * blocks_with_mulch_gt0 / total_blocks,
perc_nonNA = 100 * blocks_with_mulch_nonNA / total_blocks
)
print(block_summary)## # A tibble: 6 × 6
## IQR_Season total_blocks blocks_with_mulch_gt0 blocks_with_mulch_nonNA perc_gt0
## <chr> <int> <int> <int> <dbl>
## 1 23A 538 0 538 0
## 2 23B 499 1 0 0.200
## 3 24A 449 323 449 71.9
## 4 24B 463 1 0 0.216
## 5 25A 428 1 0 0.234
## 6 25B 441 1 0 0.227
## # ℹ 1 more variable: perc_nonNA <dbl>library(ggplot2)
ggplot(CONV_CA_Ratio,
aes(x = IQR_Season,
y = ICR_mulch_thicknesscm_planting)) +
geom_boxplot(outlier.color = "red", fill = "skyblue") +
labs(
title = "Distribution of Mulch Thickness (cm) at Planting by Season",
y = "ICR_mulch_thicknesscm_planting",
x = "IQR_Season"
) +
theme_minimal()
library(dplyr)
library(ggplot2)
data_filt <- CONV_CA_Ratio %>%
filter(!IQR_Season %in% c("23A", "23B", "24B", "25A", "25B"))
ggplot(data_filt,
aes(x = ICR_mulch_thicknesscm_planting,
y = Y_ratio,
color = crop)) +
geom_point(alpha = 0.6) +
geom_smooth(method = "lm", se = TRUE) +
labs(
title = "Yield Ratio vs Mulch Thickness (cm) at Planting by Crop (Year 23 removed)",
x = "ICR_mulch_thicknesscm_planting",
y = "Y_ratio",
color = "Crop"
) +
theme_minimal()
data_filt2 <- CONV_CA_Ratio %>%
filter(!IQR_Season %in% c("23A", "23B", "24B", "25A", "25B"))
ggplot(data_filt2,
aes(x = ICR_mulch_thicknesscm_planting,
y = ICR_mulch_perc_planti,
color = crop)) +
geom_point(alpha = 0.6) +
geom_smooth(method = "lm", se = TRUE) +
labs(
title = "% Soil mulched vs Mulch Thickness (cm) at Planting by Crop (Year 23 removed)",
x = "ICR_mulch_thicknesscm_planting",
y = "% Soil mulched",
color = "Crop"
) +
theme_minimal()
Termites
library(dplyr)
block_summary <- CONV_CA_Ratio %>%
group_by(IQR_Season) %>%
summarise(
total_blocks = n_distinct(IQR_block),
blocks_with_termite_gt0 = n_distinct(IQR_block[ICR_Termite_incidence > 0]),
blocks_with_termite_nonNA = n_distinct(IQR_block[!is.na(ICR_Termite_incidence)]),
.groups = "drop"
) %>%
mutate(
perc_gt0 = 100 * blocks_with_termite_gt0 / total_blocks,
perc_nonNA = 100 * blocks_with_termite_nonNA / total_blocks
)
print(block_summary)## # A tibble: 6 × 6
## IQR_Season total_blocks blocks_with_termite_…¹ blocks_with_termite_…² perc_gt0
## <chr> <int> <int> <int> <dbl>
## 1 23A 538 0 538 0
## 2 23B 499 1 16 0.200
## 3 24A 449 1 449 0.223
## 4 24B 463 75 463 16.2
## 5 25A 428 2 428 0.467
## 6 25B 441 1 0 0.227
## # ℹ abbreviated names: ¹blocks_with_termite_gt0, ²blocks_with_termite_nonNA
## # ℹ 1 more variable: perc_nonNA <dbl>termite_summary <- CONV_CA_Ratio %>%
group_by(IQR_Season) %>%
summarise(
count = sum(!is.na(ICR_Termite_incidence)),
min = min(ICR_Termite_incidence, na.rm = TRUE),
q25 = quantile(ICR_Termite_incidence, 0.25, na.rm = TRUE),
median = median(ICR_Termite_incidence, na.rm = TRUE),
mean = mean(ICR_Termite_incidence, na.rm = TRUE),
q75 = quantile(ICR_Termite_incidence, 0.75, na.rm = TRUE),
max = max(ICR_Termite_incidence, na.rm = TRUE),
sd = sd(ICR_Termite_incidence, na.rm = TRUE)
) %>%
arrange(IQR_Season)
print(termite_summary)## # A tibble: 6 × 9
## IQR_Season count min q25 median mean q75 max sd
## <chr> <int> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
## 1 23A 538 0 0 0 0 0 0 0
## 2 23B 16 0 0 0 0 0 0 0
## 3 24A 449 0 0 0 0.00548 0 2.46 0.116
## 4 24B 463 0 0 0 0.162 0 1 0.369
## 5 25A 428 0 0 0 0.0444 0 15 0.750
## 6 25B 0 Inf NA NA NaN NA -Inf NAHow termites relate to Y_ratio
There is not clear relation
library(dplyr)
library(ggplot2)
# -------------------------------------------------
# 1️⃣ Clean main dataset FIRST
# -------------------------------------------------
CONV_CA_Ratio_clean <- CONV_CA_Ratio %>%
filter(
!is.na(Y_ratio),
!is.na(IQR_Season),
!is.na(IQR_block)
) %>%
mutate(
ICR_Termite_incidence = as.numeric(ICR_Termite_incidence),
ICR_Termite_incidence = ifelse(is.na(ICR_Termite_incidence), 0, ICR_Termite_incidence),
ICR_mulch_perc_planti = as.numeric(ICR_mulch_perc_planti)
)
# -------------------------------------------------
# 2️⃣ Apply 24B termite replacement by block
# -------------------------------------------------
termite_24B <- CONV_CA_Ratio_clean %>%
filter(IQR_Season == "24B") %>%
dplyr::select(IQR_block, ICR_24B = ICR_Termite_incidence)
plot_data <- CONV_CA_Ratio_clean %>%
left_join(termite_24B, by = "IQR_block") %>%
mutate(
ICR_Termite_incidence = ifelse(
IQR_Season == "24B",
ICR_Termite_incidence,
ICR_24B
),
# 🔴 CRITICAL FIX
ICR_Termite_incidence = coalesce(ICR_Termite_incidence, 0),
Termite_group = ifelse(ICR_Termite_incidence > 0,
"Termite > 0",
"Termite = 0"),
Termite_group = factor(Termite_group,
levels = c("Termite = 0", "Termite > 0"))
) %>%
dplyr::select(IQR_block, IQR_Season, Y_ratio,
ICR_Termite_incidence, Termite_group, ICR_mulch_perc_planti)
# -------------------------------------------------
# 3️⃣ Summary statistics
# -------------------------------------------------
summary_stats <- plot_data %>%
group_by(IQR_Season, Termite_group) %>%
summarise(
mean_Y = mean(Y_ratio),
N = n(),
.groups = "drop"
)
# -------------------------------------------------
# 4️⃣ T-tests per season
# -------------------------------------------------
t_test_results <- plot_data %>%
group_by(IQR_Season) %>%
summarise(
p_value = tryCatch(
t.test(Y_ratio ~ Termite_group)$p.value,
error = function(e) NA_real_
),
y_pos = max(Y_ratio) * 1.05,
.groups = "drop"
)
# -------------------------------------------------
# 5️⃣ Plot
# -------------------------------------------------
ggplot(plot_data,
aes(x = Termite_group,
y = Y_ratio,
fill = Termite_group)) +
geom_boxplot(alpha = 0.7, outlier.shape = NA) +
geom_jitter(width = 0.15, alpha = 0.3) +
geom_text(data = summary_stats,
aes(y = mean_Y,
label = paste0("Mean=", round(mean_Y,2),
"\nN=", N)),
vjust = -1,
size = 3) +
geom_text(data = t_test_results,
aes(x = 1.5,
y = y_pos,
label = paste0("p = ", signif(p_value, 3))),
inherit.aes = FALSE,
size = 4) +
facet_wrap(~ IQR_Season) +
theme_minimal(base_size = 14) +
theme(
legend.position = "none",
strip.text = element_text(face = "bold")
) +
labs(
title = "Y_ratio by Termite Incidence Across Seasons",
x = "Termite Incidence",
y = "Y_ratio"
)
RElation of mulch and termites
Consistently, plots with termite have more mulch. Then while mulch related positively with CA performance, this may mask the negative effect of termites
library(dplyr)
library(ggplot2)
# -------------------------------------------------
# 1️⃣ Clean dataset
# -------------------------------------------------
CONV_CA_Ratio_clean <- CONV_CA_Ratio %>%
filter(
!is.na(ICR_mulch_perc_planti),
!is.na(IQR_Season),
!is.na(IQR_block)
) %>%
mutate(
ICR_Termite_incidence = as.numeric(ICR_Termite_incidence),
ICR_Termite_incidence = ifelse(is.na(ICR_Termite_incidence), 0, ICR_Termite_incidence),
ICR_mulch_perc_planti = as.numeric(ICR_mulch_perc_planti)
)
# -------------------------------------------------
# 2️⃣ Apply 24B termite replacement
# -------------------------------------------------
termite_24B <- CONV_CA_Ratio_clean %>%
filter(IQR_Season == "24B") %>%
dplyr::select(IQR_block, ICR_24B = ICR_Termite_incidence)
plot_data <- CONV_CA_Ratio_clean %>%
left_join(termite_24B, by = "IQR_block") %>%
mutate(
ICR_Termite_incidence = ifelse(
IQR_Season == "24B",
ICR_Termite_incidence,
ICR_24B
),
ICR_Termite_incidence = coalesce(ICR_Termite_incidence, 0),
Termite_group = ifelse(ICR_Termite_incidence > 0,
"Termite > 0",
"Termite = 0"),
Termite_group = factor(Termite_group,
levels = c("Termite = 0", "Termite > 0"))
) %>%
dplyr::select(IQR_block, IQR_Season,
ICR_mulch_perc_planti,
Termite_group)
# -------------------------------------------------
# 3️⃣ Summary statistics (Mulch)
# -------------------------------------------------
summary_stats <- plot_data %>%
group_by(IQR_Season, Termite_group) %>%
summarise(
mean_mulch = mean(ICR_mulch_perc_planti),
N = n(),
.groups = "drop"
)
# -------------------------------------------------
# 4️⃣ T-tests per season (Mulch)
# -------------------------------------------------
t_test_results <- plot_data %>%
group_by(IQR_Season) %>%
summarise(
p_value = tryCatch(
t.test(ICR_mulch_perc_planti ~ Termite_group)$p.value,
error = function(e) NA_real_
),
y_pos = max(ICR_mulch_perc_planti) * 1.05,
.groups = "drop"
)
# -------------------------------------------------
# 5️⃣ Plot
# -------------------------------------------------
ggplot(plot_data,
aes(x = Termite_group,
y = ICR_mulch_perc_planti,
fill = Termite_group)) +
geom_boxplot(alpha = 0.7, outlier.shape = NA) +
geom_jitter(width = 0.15, alpha = 0.3) +
geom_text(data = summary_stats,
aes(y = mean_mulch,
label = paste0("Mean=", round(mean_mulch,1),
"\nN=", N)),
vjust = -1,
size = 3) +
geom_text(data = t_test_results,
aes(x = 1.5,
y = y_pos,
label = paste0("p = ", signif(p_value, 3))),
inherit.aes = FALSE,
size = 4) +
facet_wrap(~ IQR_Season) +
theme_minimal(base_size = 14) +
theme(
legend.position = "none",
strip.text = element_text(face = "bold")
) +
labs(
title = "Mulch Percentage at Planting by Termite Incidence",
x = "Termite Incidence",
y = "Mulch at Planting (%)"
)
## Dive into the problematic weed “couchgrass; Digitaria abyssinica”
This weed is very competitive and hard to erradicate by hand, so very
complicated in CA systems. Current weed variable
“Weed_species_A_combined” only divides into broad leaf, grasses and
whitch weed. So now we will digg more into couch grass specifically
First lets explore when was this variable recorded NOtes: not well recorded in 23A
Extract the precence of couch grass in each block from ca_raw
library(dplyr)
library(stringr)
ca_raw <- read_csv("D:/Mes Donnees/OAF_CIRAD/CA_MoU/data/raw/ca_raw.csv")
library(dplyr)
result <- ca_raw %>%
# Keep only Treat_code C
filter(Treat_code == "C") %>%
# Detect Couch grass (any spelling variant)
mutate(Couch_Presence = ifelse(
str_detect(IQR_Weedtype1_specie,
regex("Couch_grass|couchgrass|Couch grass", ignore_case = FALSE)),
1, 0
)) %>%
# Group by block and season
group_by(IQR_block, IQR_Season) %>%
# If it appears at least once → 1
summarise(Couch_Presence = max(Couch_Presence, na.rm = TRUE),
.groups = "drop")
resultlibrary(dplyr)
season_percentages <- result %>%
group_by(IQR_Season) %>%
summarise(
pct_0 = mean(Couch_Presence == 0, na.rm = TRUE) * 100,
pct_1 = mean(Couch_Presence == 1, na.rm = TRUE) * 100,
pct_negInf = mean(Couch_Presence == -Inf, na.rm = TRUE) * 100,
.groups = "drop"
)
season_percentageslibrary(dplyr)
library(tidyr)
couch_block_table <- result %>%
# Keep only relevant seasons
filter(IQR_Season %in% c("23A", "23B", "24A", "25A", "25B")) %>%
# Pivot to wide format
pivot_wider(
names_from = IQR_Season,
values_from = Couch_Presence,
names_prefix = "Couch_"
) %>%
# Replace missing with 0 (if a block had no observation that season)
mutate(across(starts_with("Couch_"), ~replace_na(., 0))) %>%
# Create consistency variable
mutate(
Couch_Consit = ifelse(
rowSums(across(starts_with("Couch_"))) >= 3,
1, 0
)
)
couch_block_tablelibrary(dplyr)
percent_ones <- couch_block_table %>%
summarise(
across(
-IQR_block,
~ mean(. == 1, na.rm = TRUE) * 100
)
)
percent_onesNow we evaluate the relation among having couch grass (at least 2 seasons to make it more robust) with the overall average Y_ratio of a
library(dplyr)
library(ggplot2)
library(ggpubr)
plot_data <- CONV_CA_Ratio %>%
mutate(Couch_Consit = factor(Couch_Consit, levels = c(0, 1)))
# per-facet positions
y_positions <- plot_data %>%
group_by(IQF_environment) %>%
summarise(
y_max = max(Y_ratio, na.rm = TRUE),
y_label = y_max * 0.8,
y_pval = y_max * 0.6,
.groups = "drop"
)
# summary labels (N, Mean, %>=1)
summary_labels <- plot_data %>%
group_by(IQF_environment, Couch_Consit) %>%
summarise(
N = n(),
Mean = mean(Y_ratio, na.rm = TRUE),
Pct_ge1 = mean(Y_ratio >= 1, na.rm = TRUE) * 100,
label = paste0(
"N = ", N,
"\nMean = ", round(Mean, 2),
"\n≥1 = ", round(Pct_ge1, 1), "%"
),
.groups = "drop"
) %>%
left_join(y_positions, by = "IQF_environment")
# p-values per environment
pvals <- plot_data %>%
group_by(IQF_environment) %>%
summarise(p = tryCatch(t.test(Y_ratio ~ Couch_Consit)$p.value, error = function(e) NA_real_),
.groups = "drop") %>%
left_join(y_positions, by = "IQF_environment") %>%
mutate(
group1 = "0",
group2 = "1",
y.position = y_pval,
p.signif = case_when(
is.na(p) ~ "NA",
p <= 0.001 ~ "***",
p <= 0.01 ~ "**",
p <= 0.05 ~ "*",
TRUE ~ "ns"
)
)
ggplot(plot_data, aes(x = Couch_Consit, y = Y_ratio, fill = Couch_Consit)) +
geom_boxplot(alpha = 0.7, outlier.shape = NA) +
geom_text(
data = summary_labels,
aes(x = Couch_Consit, y = y_label, label = label),
inherit.aes = FALSE,
size = 3.2
) +
stat_pvalue_manual(
pvals,
label = "p.signif",
xmin = "group1",
xmax = "group2",
y.position = "y.position",
tip.length = 0.01
) +
facet_wrap(~ IQF_environment, scales = "free_y") +
labs(x = "Couch Consistency", y = "Y_ratio") +
coord_cartesian(clip = "off") +
theme_bw() +
theme(
legend.position = "none",
strip.text = element_text(face = "bold"),
plot.margin = ggplot2::margin(10, 10, 30, 10)
)
Now lets see hoe couch grass related to accumulated yield
library(dplyr)
library(ggplot2)
library(ggpubr)
# Prepare data
plot_data <- drivers_data %>%
mutate(Couch_Consit = factor(Couch_Consit, levels = c(0, 1)))
# Per-facet positions
y_positions <- plot_data %>%
group_by(IQF_environment) %>%
summarise(
y_max = max(avg_Y_ratio, na.rm = TRUE),
y_label = y_max * 0.8,
y_pval = y_max * 0.6,
.groups = "drop"
)
# Summary labels (N, Mean, %>=1)
summary_labels <- plot_data %>%
group_by(IQF_environment, Couch_Consit) %>%
summarise(
N = n(),
Mean = mean(avg_Y_ratio, na.rm = TRUE),
Pct_ge1 = mean(avg_Y_ratio >= 1, na.rm = TRUE) * 100,
label = paste0(
"N = ", N,
"\nMean = ", round(Mean, 2),
"\n≥1 = ", round(Pct_ge1, 1), "%"
),
.groups = "drop"
) %>%
left_join(y_positions, by = "IQF_environment")
# p-values per environment
pvals <- plot_data %>%
group_by(IQF_environment) %>%
summarise(
p = tryCatch(t.test(avg_Y_ratio ~ Couch_Consit)$p.value,
error = function(e) NA_real_),
.groups = "drop"
) %>%
left_join(y_positions, by = "IQF_environment") %>%
mutate(
group1 = "0",
group2 = "1",
y.position = y_pval,
p.signif = case_when(
is.na(p) ~ "NA",
p <= 0.001 ~ "***",
p <= 0.01 ~ "**",
p <= 0.05 ~ "*",
TRUE ~ "ns"
)
)
# Final plot
ggplot(plot_data, aes(x = Couch_Consit, y = avg_Y_ratio, fill = Couch_Consit)) +
geom_boxplot(alpha = 0.7, outlier.shape = NA) +
geom_text(
data = summary_labels,
aes(x = Couch_Consit, y = y_label, label = label),
inherit.aes = FALSE,
size = 3.2
) +
stat_pvalue_manual(
pvals,
label = "p.signif",
xmin = "group1",
xmax = "group2",
y.position = "y.position",
tip.length = 0.01
) +
facet_wrap(~ IQF_environment, scales = "free_y") +
labs(x = "Couch Consistency", y = "avg_Y_ratio") +
coord_cartesian(clip = "off") +
theme_bw() +
theme(
legend.position = "none",
strip.text = element_text(face = "bold"),
plot.margin = ggplot2::margin(10, 10, 30, 10)
)
###Lets see how mulch interact with couch grass precence
library(dplyr)
library(ggplot2)
# 1️⃣ Filter data
data_filt <- CONV_CA_Ratio %>%
filter(year != 23)
# 2️⃣ Fit interaction model (proper slope comparison)
model <- lm(Y_ratio ~ ICR_mulch_perc_planti * factor(Couch_Consit),
data = data_filt)
model_summary <- summary(model)
# 3️⃣ Extract slopes
coef_model <- coef(model_summary)
slope_0 <- coef_model["ICR_mulch_perc_planti", "Estimate"]
slope_diff <- coef_model["ICR_mulch_perc_planti:factor(Couch_Consit)1", "Estimate"]
slope_1 <- slope_0 + slope_diff
p_interaction <- coef_model[
"ICR_mulch_perc_planti:factor(Couch_Consit)1",
"Pr(>|t|)"
]
# 4️⃣ Create plot
ggplot(data_filt,
aes(x = ICR_mulch_perc_planti,
y = Y_ratio,
color = factor(Couch_Consit))) +
geom_point(alpha = 0.6) +
geom_smooth(method = "lm", se = TRUE) +
annotate("text",
x = Inf,
y = Inf,
hjust = 1.1,
vjust = 1.2,
label = paste0(
"Slope (0) = ", round(slope_0, 3),
"\nSlope (1) = ", round(slope_1, 3),
"\nSlope difference p = ",
signif(p_interaction, 3)
)) +
labs(
title = "Yield Ratio vs Mulch Percentage per Plant\n(Year 23 removed)",
x = "ICR_mulch_perc_planti",
y = "Y_ratio",
color = "Couch_Consit"
) +
theme_minimal()
## Now we evaluate better all the rest of the weeds
library(dplyr)
library(tidyr)
weeds_data <- CONV_CA_Ratio %>%
dplyr::select(IQR_block, IQR_Season, IQF_Weed_specie) %>%
separate(
IQF_Weed_specie,
into = paste0("Weed_", 1:13
),
sep = " ",
fill = "right",
extra = "drop"
)Now will see the list of unique weeds to correct if there is errors
library(dplyr)
library(tidyr)
unique_weeds <- weeds_data %>%
pivot_longer(
cols = starts_with("Weed_"),
names_to = "Weed_number",
values_to = "Weed_name"
) %>%
filter(!is.na(Weed_name)) %>%
distinct(Weed_name) %>%
arrange(Weed_name)
unique_weedsNow we add a variable that counts the N of weeds species in the plot
library(dplyr)
weeds_data <- weeds_data %>%
rowwise() %>%
mutate(
N_weeds = n_distinct(
c_across(starts_with("Weed_"))[
!c_across(starts_with("Weed_")) %in% c(NA, "Other", "--")
]
)
) %>%
ungroup()Lets see how was the weeds N across season
library(ggplot2)
library(dplyr)
ggplot(weeds_data, aes(x = N_weeds)) +
geom_bar(fill = "steelblue", color = "black") +
facet_wrap(~ IQR_Season) +
theme_bw() +
labs(
x = "Number of Unique Weeds (N_weeds)",
y = "Frequency (Number of Plots)",
title = "Frequency Distribution of Weed Richness by Season"
)
Now we make a new variable that indicats weed pressure
1 if 1 N_weeds = 0 1.5 if N_Broadleaf = 1 & N_Grass = 0 2 if N_Broadleaf = >2 & N_Grass = 0 3 if N_Broadleaf = 0 & N_Grass = 1 & CouchGrass is = 0 3.5 if N_Broadleaf > 0 & N_Grass = 1 & CouchGrass is = 0 4 N_Grass > 1 & CouchGrass is = 0 6 if N_Broadleaf = 0 & N_Grass = 1 & CouchGrass is = 1 6.5 if N_Broadleaf > 0 & N_Grass = 1 & CouchGrass is = 1 8 if N_Grass > 1 & CouchGrass is = 1 ### Same but for which weed ### now lets add the variable in drivers_data ### Lets akso add the variable Whitch_weed into drivers data
lets see how whitch weed related to CA performance
library(ggplot2)
library(dplyr)
# Convert to factor once, up front, so everything downstream is consistent
drivers_data <- drivers_data %>%
mutate(Whitch_weed = factor(Whitch_weed, levels = c(0, 1)))
# Create summary stats
summary_stats <- drivers_data %>%
group_by(IQF_environment, Whitch_weed) %>%
summarise(
N = n(),
Mean = mean(avg_Y_ratio, na.rm = TRUE),
y_pos = max(avg_Y_ratio, na.rm = TRUE) * 1.05,
.groups = "drop"
)
# Plot
ggplot(drivers_data, aes(x = Whitch_weed, y = avg_Y_ratio, fill = Whitch_weed)) +
geom_boxplot(alpha = 0.7) +
# Mean point inside box
stat_summary(fun = mean, geom = "point",
shape = 23, size = 3, fill = "white") +
# Add N and Mean above box
geom_text(
data = summary_stats,
aes(x = Whitch_weed, y = y_pos,
label = paste0("N=", N, "\nMean=", round(Mean, 2))),
size = 3,
inherit.aes = FALSE
) +
facet_wrap(~ IQF_environment) +
scale_fill_manual(values = c("0" = "steelblue", "1" = "tomato")) +
theme_bw() +
theme(
legend.position = "none",
strip.text = element_text(face = "bold"),
plot.title = element_text(hjust = 0.5)
) +
labs(
x = "Whitch Weed Presence",
y = "Y Ratio",
title = "Y_ratio by Whitch Weed Presence Across Environments"
)
###
library(ggplot2)
library(dplyr)
library(ggpubr)
# Summary stats
summary_stats <- drivers_data %>%
group_by(Whitch_weed) %>%
summarise(
N = n(),
Mean = mean(avg_Y_ratio, na.rm = TRUE),
y_pos = max(avg_Y_ratio, na.rm = TRUE) * 1.08,
.groups = "drop"
)
# Plot
ggplot(drivers_data, aes(x = Whitch_weed, y = avg_Y_ratio, fill = Whitch_weed)) +
geom_boxplot(alpha = 0.7) +
stat_summary(fun = mean, geom = "point",
shape = 23, size = 3, fill = "white") +
# FIX: include x aesthetic
geom_text(
data = summary_stats,
aes(x = Whitch_weed,
y = y_pos,
label = paste0("N=", N, "\nMean=", round(Mean, 2))),
size = 4
) +
stat_compare_means(method = "t.test",
label = "p.format") +
scale_fill_manual(values = c("0" = "steelblue",
"1" = "tomato")) +
theme_bw() +
theme(
legend.position = "none",
plot.title = element_text(hjust = 0.5)
) +
labs(
x = "Whitch Weed Presence",
y = "Average Y Ratio",
title = "Effect of Whitch Weed Presence on Y_ratio"
)
library(dplyr)
percent_blocks <- drivers_data %>%
# Reduce to unique block-environment combinations
distinct(IQF_environment, IQR_block, avg_Y_ratio) %>%
# Create indicator per block
mutate(ge1 = avg_Y_ratio >= 1) %>%
# Summarise per environment
group_by(IQF_environment) %>%
summarise(
n_blocks = n(),
n_ge1 = sum(ge1, na.rm = TRUE),
percent_ge1 = 100 * n_ge1 / n_blocks,
.groups = "drop"
)
percent_blockslibrary(dplyr)
overall_percent <- drivers_data %>%
# Ensure one row per block
distinct(IQR_block, avg_Y_ratio) %>%
summarise(
n_blocks = n(),
n_ge1 = sum(avg_Y_ratio>= 1, na.rm = TRUE),
percent_ge1 = 100 * n_ge1 / n_blocks
)
overall_percentSH residues at harvest
library(dplyr)
block_summary <- CONV_CA_Ratio %>%
group_by(IQR_Season) %>%
summarise(
total_blocks = n_distinct(IQR_block),
blocks_with_sh_residues_gt0 = n_distinct(IQR_block[ICR_sh_residues_tn_ha > 0]),
blocks_with_sh_residues_nonNA = n_distinct(IQR_block[!is.na(ICR_sh_residues_tn_ha)]),
.groups = "drop"
) %>%
mutate(
perc_gt0 = 100 * blocks_with_sh_residues_gt0 / total_blocks,
perc_nonNA = 100 * blocks_with_sh_residues_nonNA / total_blocks
)
print(block_summary)## # A tibble: 6 × 6
## IQR_Season total_blocks blocks_with_sh_resid…¹ blocks_with_sh_resid…² perc_gt0
## <chr> <int> <int> <int> <dbl>
## 1 23A 538 0 538 0
## 2 23B 499 1 0 0.200
## 3 24A 449 303 449 67.5
## 4 24B 463 199 463 43.0
## 5 25A 428 278 428 65.0
## 6 25B 441 1 0 0.227
## # ℹ abbreviated names: ¹blocks_with_sh_residues_gt0,
## # ²blocks_with_sh_residues_nonNA
## # ℹ 1 more variable: perc_nonNA <dbl>sh_residue_summary <- CONV_CA_Ratio %>%
group_by(IQR_Season) %>%
summarise(
count = sum(!is.na(ICR_sh_residues_tn_ha)),
min = min(ICR_sh_residues_tn_ha, na.rm = TRUE),
q25 = quantile(ICR_sh_residues_tn_ha, 0.25, na.rm = TRUE),
median = median(ICR_sh_residues_tn_ha, na.rm = TRUE),
mean = mean(ICR_sh_residues_tn_ha, na.rm = TRUE),
q75 = quantile(ICR_sh_residues_tn_ha, 0.75, na.rm = TRUE),
max = max(ICR_sh_residues_tn_ha, na.rm = TRUE),
sd = sd(ICR_sh_residues_tn_ha, na.rm = TRUE)
) %>%
arrange(IQR_Season)
print(sh_residue_summary)## # A tibble: 6 × 9
## IQR_Season count min q25 median mean q75 max sd
## <chr> <int> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
## 1 23A 538 0 0 0 0 0 0 0
## 2 23B 0 Inf NA NA NaN NA -Inf NA
## 3 24A 449 0 0 3.21 4.57 7.23 33.1 5.21
## 4 24B 463 0 0 0 0.785 0.695 12 1.63
## 5 25A 428 0 0 1.01 3.96 4.61 120 9.07
## 6 25B 0 Inf NA NA NaN NA -Inf NAlibrary(ggplot2)
ggplot(CONV_CA_Ratio,
aes(x = IQR_Season,
y = ICR_sh_residues_tn_ha)) +
geom_boxplot(outlier.color = "red", fill = "skyblue") +
labs(
title = "Distribution of SH Residues (t/ha) by Season",
y = "ICR_sh_residues_tn_ha",
x = "IQR_Season"
) +
theme_minimal()
model <- lm(Y_ratio ~ ICR_sh_residues_tn_ha,
data = CONV_CA_Ratio,
na.action = na.omit)
summary(model)##
## Call:
## lm(formula = Y_ratio ~ ICR_sh_residues_tn_ha, data = CONV_CA_Ratio,
## na.action = na.omit)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.73520 -0.19663 -0.02359 0.15158 2.93672
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 0.8013764 0.0079983 100.193 <2e-16 ***
## ICR_sh_residues_tn_ha -0.0002543 0.0013606 -0.187 0.852
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.3217 on 1876 degrees of freedom
## (940 observations deleted due to missingness)
## Multiple R-squared: 1.862e-05, Adjusted R-squared: -0.0005144
## F-statistic: 0.03494 on 1 and 1876 DF, p-value: 0.8518library(dplyr)
library(ggplot2)
data_filt <- CONV_CA_Ratio %>%
filter(!IQR_Season %in% c("23A", "23B", "24B", "25B"))
ggplot(data_filt,
aes(x = ICR_sh_residues_tn_ha,
y = Y_ratio,
color = crop)) +
geom_point(alpha = 0.6) +
geom_smooth(method = "lm", se = TRUE) +
labs(
title = "Yield Ratio vs SH Residues (t/ha) by Crop (Year 23 removed)",
x = "ICR_sh_residues_tn_ha",
y = "Y_ratio",
color = "Crop"
) +
theme_minimal()