R code for At-plant fertilizer strategy and soil organic N cycling govern early-season maize N uptake and fertilizer N retention
Manjot Kaur Rekhi
2026-10-25
Load libraries
library(car)## Loading required package: carDatarequire(ggplot2)## Loading required package: ggplot2require(tidyverse)## Loading required package: tidyverse## ── Attaching core tidyverse packages ──────────────────────── tidyverse 2.0.0 ──
## ✔ dplyr 1.2.1 ✔ readr 2.2.0
## ✔ forcats 1.0.1 ✔ stringr 1.6.0
## ✔ lubridate 1.9.5 ✔ tibble 3.3.1
## ✔ purrr 1.2.2 ✔ tidyr 1.3.2
## ── Conflicts ────────────────────────────────────────── tidyverse_conflicts() ──
## ✖ dplyr::filter() masks stats::filter()
## ✖ dplyr::lag() masks stats::lag()
## ✖ dplyr::recode() masks car::recode()
## ✖ purrr::some() masks car::some()
## ℹ Use the conflicted package (<http://conflicted.r-lib.org/>) to force all conflicts to become errorslibrary(tidyr)
library(dplyr)
library(emmeans)## Welcome to emmeans.
## Caution: You lose important information if you filter this package's results.
## See '? untidy'library(agricolae)
library(multcomp)## Loading required package: mvtnorm
## Loading required package: survival
## Loading required package: TH.data
## Loading required package: MASS
##
## Attaching package: 'MASS'
##
## The following object is masked from 'package:dplyr':
##
## select
##
##
## Attaching package: 'TH.data'
##
## The following object is masked from 'package:MASS':
##
## geyserlibrary(multcompView)
library(ggpubr)
library(ggpmisc)## Loading required package: ggpp
## Registered S3 methods overwritten by 'ggpp':
## method from
## heightDetails.titleGrob ggplot2
## widthDetails.titleGrob ggplot2
##
## Attaching package: 'ggpp'
##
## The following objects are masked from 'package:ggpubr':
##
## as_npc, as_npcx, as_npcy
##
## The following object is masked from 'package:ggplot2':
##
## annotatelibrary(ggthemes)
library(nlme)##
## Attaching package: 'nlme'
##
## The following object is masked from 'package:dplyr':
##
## collapselibrary(lme4)## Loading required package: Matrix
##
## Attaching package: 'Matrix'
##
## The following objects are masked from 'package:tidyr':
##
## expand, pack, unpack
##
## Registered S3 method overwritten by 'lme4':
## method from
## na.action.merMod car
##
## Attaching package: 'lme4'
##
## The following object is masked from 'package:nlme':
##
## lmListlibrary(lmerTest)##
## Attaching package: 'lmerTest'
##
## The following object is masked from 'package:lme4':
##
## lmer
##
## The following object is masked from 'package:stats':
##
## steplibrary(readxl)
library(openxlsx)Import Data
#import data
npj_data<-read_excel("npj_data.xlsx")
#convert variables to factors
npj_data$location<-as.factor(npj_data$location)
npj_data$row<-as.factor(npj_data$row)
npj_data$trt<-as.factor(npj_data$trt)
npj_data$position<-as.factor(npj_data$position)
npj_data$stage<-as.factor(npj_data$stage)
str(npj_data)## tibble [192 × 45] (S3: tbl_df/tbl/data.frame)
## $ stage : Factor w/ 2 levels "V3","V6": 1 1 1 1 1 1 1 1 1 1 ...
## $ location : Factor w/ 2 levels "Neb","NH": 1 1 1 1 1 1 1 1 1 1 ...
## $ row : Factor w/ 8 levels "R1","R2","R3",..: 1 1 1 1 1 1 3 3 3 3 ...
## $ trt : Factor w/ 6 levels "Both 1x","Both 2x",..: 3 4 5 6 1 2 3 4 5 6 ...
## $ position : Factor w/ 2 levels "Band","In-row": 2 2 2 2 2 2 2 2 2 2 ...
## $ totN_applied_kgha : num [1:192] 5.18 10.36 73.98 147.95 79.15 ...
## $ avg_shoot_wt_g : num [1:192] 0.823 0.77 0.463 0.395 0.39 ...
## $ avg_root_wt_g : num [1:192] 0.4 0.25 0.26 0.2 0.26 0.39 0.36 0.34 0.2 0.28 ...
## $ avg_root_wt_band : num [1:192] NA NA NA NA NA NA NA NA NA NA ...
## $ avg_root_wt_tot : num [1:192] 0.4 0.25 0.26 0.2 0.26 0.39 0.36 0.34 0.2 0.28 ...
## $ totN_shoot_kgha : num [1:192] 2.95 2.9 2.11 2 1.69 ...
## $ totN_root_kgha : num [1:192] 0.705 0.504 0.595 0.539 0.593 ...
## $ APE_shoot : num [1:192] 2.44 5.43 2.45 3.99 3.18 ...
## $ APE_root : num [1:192] 2.66 3.83 2.02 2.45 2.96 ...
## $ perNdff_shoot : num [1:192] 29.2 65 29.4 47.8 38.1 ...
## $ perN_dff_root : num [1:192] 31.9 45.9 24.2 29.3 35.5 ...
## $ FertN_shoot_kg_ha : num [1:192] 0.863 1.888 0.619 0.956 0.643 ...
## $ FertN_root_kg_ha : num [1:192] 0.225 0.231 0.144 0.158 0.211 ...
## $ perRecov_shoot : num [1:192] 16.671 18.23 0.836 0.646 0.812 ...
## $ perRecov_root : num [1:192] 4.335 2.227 0.194 0.106 0.265 ...
## $ totC_shoot_kgha : num [1:192] 30.8 29.2 17.2 14.6 14.9 ...
## $ totC_root_kgha : num [1:192] 14.95 9.21 9.3 7.33 9.59 ...
## $ ug_NO3-N_g : num [1:192] 0.753 2.583 1.526 2.257 1.62 ...
## $ ug_NH4-N_g : num [1:192] 3.729 6.263 1.243 0.895 4.393 ...
## $ ug_DIN_g : num [1:192] 4.48 8.85 2.77 3.15 6.01 ...
## $ pH : num [1:192] 5.24 6.44 6.3 6.17 5.94 6.08 6.02 6.5 6.58 6.44 ...
## $ fracMAOM_gperg : num [1:192] 0.859 0.879 0.85 0.835 0.849 ...
## $ fracPOM_gperg : num [1:192] 0.061 0.041 0.0972 0.082 0.0699 ...
## $ perNdff_TDN : num [1:192] 0.459 12.593 0.572 0.678 1.739 ...
## $ kgN_dff_TDN_ha : num [1:192] 0.00322 0.11166 0.00349 0.0045 0.01102 ...
## $ kgN_dff_TDN_ha_band : num [1:192] 0.01573 0.00387 1.39498 4.55123 8.20886 ...
## $ kgN_dff_TDN_ha_furrowband : num [1:192] 0.0189 0.1155 1.3985 4.5557 8.2199 ...
## $ per_Recov_TDN_tot : num [1:192] 0.366 1.115 1.89 3.079 10.385 ...
## $ perN_dff_wholesoil : num [1:192] 0.502 2.091 0.158 0.69 1.366 ...
## $ kgN_dff_wholesoil_ha : num [1:192] 0.888 3.646 0.282 1.294 2.479 ...
## $ kgN_dff_wholesoil_ha_furrowband: num [1:192] 1.11 4.05 15.55 48.34 23.81 ...
## $ per_Recov_wholesoil_tot : num [1:192] 21.5 39.1 21 32.7 30.1 ...
## $ perN_dff_MAOM : num [1:192] 0.372 0.985 0.136 0.43 0.72 ...
## $ kgN_dff_MAOM_ha : num [1:192] 0.547 1.427 0.206 0.595 0.95 ...
## $ kgN_dff_MAOM_ha_furrowband : num [1:192] 0.83 1.86 8.29 17.69 9.88 ...
## $ per_Recov_MAOM_tot : num [1:192] 16 18 11.2 12 12.5 ...
## $ kgN_dff_MBN_ha : num [1:192] 0.01473 0 0.00101 0.0237 0.00318 ...
## $ kgN_dff_MBN_ha_band : num [1:192] 0 0.00558 0.44992 0.54043 1.07801 ...
## $ kgN_dff_MBN_ha_furrowband : num [1:192] 0.01473 0.00558 0.45093 0.56414 1.08119 ...
## $ per_Recov_MBN_tot : num [1:192] 0.2844 0.0539 0.6096 0.3813 1.3659 ...Data wrangling
# Duplicate data
data_new <- npj_data
# Reorder factor levels
data_new$trt <- factor(data_new$trt,
c( "InF 1x", "InF 2x", "UAN 1x"
, "Both 1x","UAN 2x", "Both 2x"))
#Recode trt names
data_new<- data_new %>%
mutate(trt = dplyr::recode(
trt,
"InF 1x" = "InF-5",
"InF 2x" = "InF-10",
"UAN 1x" = "UAN-74",
"Both 1x" = "Both-79",
"UAN 2x" = "UAN-148",
"Both 2x" = "Both-158"
))
#Reorder factor levels
data_new$trt <- factor(data_new$trt,
c( "InF-5", "InF-10", "UAN-74"
, "Both-79","UAN-148", "Both-158"))
#Recode location names
data_new <- data_new |>
dplyr::mutate(location = dplyr::recode(location,
"Neb" = "NE",
"NH" = "NH"
))
#Convert root and shoot weight from g to kg/ha and make new columns
data_new$shoot_wt<- (data_new$avg_shoot_wt_g/1000)*86000
data_new$root_wt<- (data_new$avg_root_wt_g/1000)*86000
data_new$root_wt_tot<- (data_new$avg_root_wt_tot/1000)*86000
#only furrow roots
data_new$root_wt_inf <- ifelse(
data_new$position == "In-row",
(data_new$avg_root_wt_g / 1000) * 86000,
NA_real_
)
data_new$totN_root_kgha_inf <- ifelse(
data_new$position == "In-row", data_new$totN_root_kgha, NA_real_
)
data_new$FertN_root_kg_ha_inf <- ifelse(
data_new$position == "In-row", data_new$FertN_root_kg_ha, NA_real_
)
data_new$perRecov_root_inf <- ifelse(
data_new$position == "In-row", data_new$perRecov_root, NA_real_
)
#Calculate soil-derived N (kg/ha) in shoot and root and make new columns
data_new$Ndf_notfert<-(data_new$totN_shoot_kgha)-(data_new$FertN_shoot_kg_ha)
data_new$Ndf_notfertroot<-(data_new$totN_root_kgha)-(data_new$FertN_root_kg_ha)
#SUBSET DATA
dat_Neb<-subset(data_new, data_new$location == "NE")
dat_NH<-subset(data_new, data_new$location == "NH")
dat_V3<-subset(data_new, data_new$stage == "V3")
dat_V6<-subset(data_new, data_new$stage == "V6")
dat_Neb_V3<-subset(dat_Neb, dat_Neb$stage == "V3")
dat_NH_V3<-subset(dat_NH, dat_NH$stage == "V3")
dat_Neb_V6<-subset(dat_Neb, dat_Neb$stage == "V6")
dat_NH_V6<-subset(dat_NH, dat_NH$stage == "V6")
dat_Neb_V3_Inf<-subset(dat_Neb_V3, dat_Neb_V3$position == "In-row")
dat_NH_V3_Band<-subset(dat_NH_V3, dat_NH_V3$position== "Band")
dat_Neb_V3_Band<-subset(dat_Neb_V3, dat_Neb_V3$position == "Band")
dat_NH_V3_Inf<-subset(dat_NH_V3, dat_NH_V3$position== "In-row")
dat_Neb_V6_Inf<-subset(dat_Neb_V6, dat_Neb_V6$position == "In-row")
dat_NH_V6_Band<-subset(dat_NH_V6, dat_NH_V6$position== "Band")
dat_Neb_V6_Band<-subset(dat_Neb_V6, dat_Neb_V6$position == "Band")
dat_NH_V6_Inf<-subset(dat_NH_V6, dat_NH_V6$position== "In-row")
InF<-subset(data_new, data_new$position == "In-row")
Band<-subset(data_new, data_new$position == "Band")
InF_V3<-subset(InF, InF$stage == "V3")
Band_V3<-subset(Band, Band$stage == "V3")
InF_V6<-subset(InF, InF$stage == "V6")
Band_V6<-subset(Band, Band$stage == "V6")Theme Manjot
theme_Manjot <- function(base_size=14, base_family="TT Times New Roman") {(theme_foundation(base_size=base_size, base_family=base_family)+
theme_bw()+
theme(
panel.background = element_rect(colour = NA),
plot.background = element_rect(colour = NA),
axis.title = element_text(color="black",size=rel(1.5)),
axis.text = element_text(color="black", size =rel(1.2)),
legend.key = element_rect(colour = NA),
legend.spacing = unit(0, "cm"),
legend.text = element_text(size= rel(1.4)),
legend.title = element_text(size=rel(1.4)),
panel.grid = element_blank(),
plot.title = element_text(color="Black",size = rel(2),face = "bold",hjust = 0),
strip.text = element_text(color="Black",size = rel(1.5),face="bold")
))
}Overall ANOVA
#variables with both position data separate
soil_vars <- c(
"kgN_dff_MAOM_ha",
"kgN_dff_wholesoil_ha",
"kgN_dff_TDN_ha",
"kgN_dff_MBN_ha",
"ug_NO3-N_g",
"ug_NH4-N_g",
"ug_DIN_g",
"pH"
)
#variable with combined/single position data
plant_vars <- c(
"shoot_wt",
"root_wt_tot",
"root_wt_inf",
"totN_shoot_kgha",
"FertN_shoot_kg_ha",
"totN_root_kgha_inf",
"FertN_root_kg_ha_inf",
"perRecov_shoot",
"perRecov_root_inf",
"Ndf_notfert",
"per_Recov_TDN_tot",
"per_Recov_MBN_tot",
"per_Recov_wholesoil_tot",
"per_Recov_MAOM_tot"
)
get_anova_soil <- function(response_var, data) {
dat <- data %>%
dplyr::select(all_of(c(response_var, "trt", "location", "stage", "position", "row"))) %>%
dplyr::filter(!is.na(.data[[response_var]])) %>%
droplevels()
if (
nrow(dat) == 0 ||
nlevels(factor(dat$trt)) < 2 ||
nlevels(factor(dat$location)) < 2 ||
nlevels(factor(dat$stage)) < 2 ||
nlevels(factor(dat$position)) < 2
) {
return(
tibble(
response = response_var,
term = NA_character_,
p_value = NA_real_,
model_type = "soil_full",
status = "skipped"
)
)
}
form <- as.formula(
paste(response_var, "~ trt * location * stage * position + (1 | location:row)")
)
mod <- tryCatch(
lmer(form, data = dat),
error = function(e) e
)
if (inherits(mod, "error")) {
return(
tibble(
response = response_var,
term = NA_character_,
p_value = NA_real_,
model_type = "soil_full",
status = mod$message
)
)
}
anova(mod) %>%
as.data.frame() %>%
tibble::rownames_to_column("term") %>%
dplyr::transmute(
response = response_var,
term = term,
p_value = `Pr(>F)`,
model_type = "soil_full",
status = "ok"
)
}
get_anova_plant <- function(response_var, data) {
dat <- data %>%
dplyr::select(all_of(c(response_var, "trt", "location", "stage", "row"))) %>%
dplyr::filter(!is.na(.data[[response_var]])) %>%
droplevels()
if (
nrow(dat) == 0 ||
nlevels(factor(dat$trt)) < 2 ||
nlevels(factor(dat$location)) < 2 ||
nlevels(factor(dat$stage)) < 2
) {
return(
tibble(
response = response_var,
term = NA_character_,
p_value = NA_real_,
model_type = "plant_no_position",
status = "skipped"
)
)
}
form <- as.formula(
paste(response_var, "~ trt * location * stage + (1 | location:row)")
)
mod <- tryCatch(
lmer(form, data = dat),
error = function(e) e
)
if (inherits(mod, "error")) {
return(
tibble(
response = response_var,
term = NA_character_,
p_value = NA_real_,
model_type = "plant_no_position",
status = mod$message
)
)
}
anova(mod) %>%
as.data.frame() %>%
tibble::rownames_to_column("term") %>%
dplyr::transmute(
response = response_var,
term = term,
p_value = `Pr(>F)`,
model_type = "plant_no_position",
status = "ok"
)
}
soil_anova_table <- purrr::map_dfr(soil_vars, get_anova_soil, data = data_new)## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')plant_anova_table <- purrr::map_dfr(plant_vars, get_anova_plant, data = data_new)## boundary (singular) fit: see help('isSingular')final_anova_table <- dplyr::bind_rows(soil_anova_table, plant_anova_table)
final_table_wide <- final_anova_table %>%
dplyr::filter(!is.na(term)) %>%
dplyr::select(response, model_type, term, p_value) %>%
tidyr::pivot_wider(
names_from = term,
values_from = p_value
)
final_table_wide <- final_table_wide %>%
dplyr::mutate(across(where(is.numeric), ~ case_when(
. < 0.001 ~ "<0.001",
TRUE ~ sprintf("%.3f", .)
)))
write.xlsx(final_table_wide, "final_anova_table14.xlsx")## Warning in file.create(to[okay]): cannot create file
## 'final_anova_table14.xlsx', reason 'Permission denied'knitr::kable(
final_table_wide,
format = "html",
caption = "ANOVA full dataset"
)| response | model_type | trt | location | stage | position | trt:location | trt:stage | location:stage | trt:position | location:position | stage:position | trt:location:stage | trt:location:position | trt:stage:position | location:stage:position | trt:location:stage:position |
|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|
| kgN_dff_MAOM_ha | soil_full | <0.001 | 0.023 | 0.001 | <0.001 | 0.867 | 0.725 | 0.133 | <0.001 | 0.285 | 0.006 | 0.697 | 0.840 | 0.136 | 0.556 | 0.827 |
| kgN_dff_wholesoil_ha | soil_full | <0.001 | <0.001 | <0.001 | <0.001 | 0.032 | 0.010 | 0.290 | <0.001 | <0.001 | <0.001 | 0.808 | 0.034 | 0.002 | 0.312 | 0.881 |
| kgN_dff_TDN_ha | soil_full | <0.001 | 0.002 | <0.001 | <0.001 | <0.001 | <0.001 | 0.118 | <0.001 | <0.001 | <0.001 | 0.266 | <0.001 | <0.001 | 0.014 | 0.234 |
| kgN_dff_MBN_ha | soil_full | 0.063 | 0.024 | 0.135 | <0.001 | 0.269 | 0.613 | 0.552 | 0.064 | 0.005 | 0.094 | 0.965 | 0.259 | 0.638 | 0.490 | 0.959 |
| ug_DIN_g | soil_full | <0.001 | <0.001 | <0.001 | <0.001 | 0.002 | <0.001 | 0.018 | <0.001 | <0.001 | <0.001 | 0.129 | 0.002 | <0.001 | 0.003 | 0.123 |
| pH | soil_full | <0.001 | <0.001 | 0.849 | 0.399 | 0.001 | <0.001 | 0.928 | 0.002 | 0.339 | <0.001 | 0.161 | 0.021 | <0.001 | <0.001 | 0.769 |
| shoot_wt | plant_no_position | 0.130 | 0.001 | <0.001 | NA | 0.376 | 0.428 | 0.002 | NA | NA | NA | 0.377 | NA | NA | NA | NA |
| root_wt_tot | plant_no_position | 0.023 | 0.007 | <0.001 | NA | 0.253 | 0.071 | 0.004 | NA | NA | NA | 0.178 | NA | NA | NA | NA |
| root_wt_inf | plant_no_position | 0.048 | 0.062 | <0.001 | NA | 0.547 | 0.292 | 0.028 | NA | NA | NA | 0.421 | NA | NA | NA | NA |
| totN_shoot_kgha | plant_no_position | <0.001 | <0.001 | <0.001 | NA | 0.116 | <0.001 | <0.001 | NA | NA | NA | 0.234 | NA | NA | NA | NA |
| FertN_shoot_kg_ha | plant_no_position | <0.001 | <0.001 | <0.001 | NA | <0.001 | <0.001 | <0.001 | NA | NA | NA | <0.001 | NA | NA | NA | NA |
| totN_root_kgha_inf | plant_no_position | <0.001 | 0.020 | <0.001 | NA | 0.103 | <0.001 | 0.004 | NA | NA | NA | 0.003 | NA | NA | NA | NA |
| FertN_root_kg_ha_inf | plant_no_position | <0.001 | 0.005 | <0.001 | NA | 0.182 | <0.001 | <0.001 | NA | NA | NA | 0.008 | NA | NA | NA | NA |
| perRecov_shoot | plant_no_position | <0.001 | 0.034 | <0.001 | NA | 0.024 | 0.236 | 0.025 | NA | NA | NA | 0.150 | NA | NA | NA | NA |
| perRecov_root_inf | plant_no_position | <0.001 | 0.003 | 0.037 | NA | <0.001 | 0.970 | 0.003 | NA | NA | NA | 0.008 | NA | NA | NA | NA |
| Ndf_notfert | plant_no_position | 0.240 | <0.001 | <0.001 | NA | 0.138 | 0.402 | <0.001 | NA | NA | NA | 0.091 | NA | NA | NA | NA |
| per_Recov_TDN_tot | plant_no_position | <0.001 | 0.001 | <0.001 | NA | <0.001 | <0.001 | 0.154 | NA | NA | NA | 0.424 | NA | NA | NA | NA |
| per_Recov_MBN_tot | plant_no_position | 0.052 | 0.039 | 0.187 | NA | 0.144 | 0.740 | 0.557 | NA | NA | NA | 0.938 | NA | NA | NA | NA |
| per_Recov_wholesoil_tot | plant_no_position | 0.007 | <0.001 | <0.001 | NA | 0.138 | 0.002 | 0.866 | NA | NA | NA | 0.477 | NA | NA | NA | NA |
| per_Recov_MAOM_tot | plant_no_position | <0.001 | 0.003 | 0.001 | NA | <0.001 | <0.001 | 0.736 | NA | NA | NA | 0.009 | NA | NA | NA | NA |
Crop biomass graph
#Shoot biomass kg ha-1 model
model_v3_neb_shoot <- lm(shoot_wt ~ trt + row, data = dat_Neb_V3)
model_v6_neb_shoot <- lm(shoot_wt ~ trt + row, data = dat_Neb_V6)
model_v3_nh_shoot <- lm(shoot_wt ~ trt + row, data = dat_NH_V3)
model_v6_nh_shoot <- lm(shoot_wt ~ trt + row, data = dat_NH_V6)
emm_v3_neb_shoot <- emmeans(model_v3_neb_shoot, ~ trt)
emm_v6_neb_shoot <- emmeans(model_v6_neb_shoot, ~ trt)
emm_v3_nh_shoot <- emmeans(model_v3_nh_shoot, ~ trt)
emm_v6_nh_shoot <- emmeans(model_v6_nh_shoot, ~ trt)
cld_v3_neb_shoot <- cld(emm_v3_neb_shoot, Letters = letters, alpha = 0.05, sort= TRUE, reversed= TRUE)
cld_v6_neb_shoot <- cld(emm_v6_neb_shoot, Letters = letters, alpha = 0.05, sort= TRUE, reversed= TRUE)
cld_v3_nh_shoot <- cld(emm_v3_nh_shoot, Letters = letters, alpha = 0.05, sort= TRUE, reversed= TRUE)
cld_v6_nh_shoot <- cld(emm_v6_nh_shoot, Letters = letters, alpha = 0.05, sort= TRUE, reversed= TRUE)
shoot_cld_df <- bind_rows(
as.data.frame(cld_v3_neb_shoot) %>% mutate(stage = "V3", location = "NE"),
as.data.frame(cld_v6_neb_shoot) %>% mutate(stage = "V6", location = "NE"),
as.data.frame(cld_v3_nh_shoot) %>% mutate(stage = "V3", location = "NH"),
as.data.frame(cld_v6_nh_shoot) %>% mutate(stage = "V6", location = "NH")
) %>%
dplyr::rename(shoot_cld = .group) %>%
dplyr::mutate(shoot_cld = gsub(" ", "", shoot_cld)) %>%
dplyr::select(trt, location, stage, shoot_cld)
# ROOT models
model_v3_neb_root <- lm(root_wt_inf ~ trt + row, data = dat_Neb_V3_Inf)
model_v6_neb_root <- lm(root_wt_inf ~ trt + row, data = dat_Neb_V6_Inf)
model_v3_nh_root <- lm(root_wt_inf ~ trt + row, data = dat_NH_V3_Inf)
model_v6_nh_root <- lm(root_wt_inf ~ trt + row, data = dat_NH_V6_Inf)
emm_v3_neb_root <- emmeans(model_v3_neb_root, ~ trt)
emm_v6_neb_root <- emmeans(model_v6_neb_root, ~ trt)
emm_v3_nh_root <- emmeans(model_v3_nh_root, ~ trt)
emm_v6_nh_root <- emmeans(model_v6_nh_root, ~ trt)
cld_v3_neb_root <- cld(emm_v3_neb_root, Letters = letters, alpha = 0.05, sort= TRUE, reversed= TRUE)
cld_v6_neb_root <- cld(emm_v6_neb_root, Letters = letters, alpha = 0.05, sort= TRUE, reversed= TRUE)
cld_v3_nh_root <- cld(emm_v3_nh_root, Letters = letters, alpha = 0.05, sort= TRUE, reversed= TRUE)
cld_v6_nh_root <- cld(emm_v6_nh_root, Letters = letters, alpha = 0.05, sort= TRUE, reversed= TRUE)
root_cld_df <- bind_rows(
as.data.frame(cld_v3_neb_root) %>% mutate(stage = "V3", location = "NE"),
as.data.frame(cld_v6_neb_root) %>% mutate(stage = "V6", location = "NE"),
as.data.frame(cld_v3_nh_root) %>% mutate(stage = "V3", location = "NH"),
as.data.frame(cld_v6_nh_root) %>% mutate(stage = "V6", location = "NH")
) %>%
dplyr:: rename(root_cld = .group) %>%
dplyr:: mutate(root_cld = gsub(" ", "", root_cld)) %>%
dplyr:: select(trt, location, stage, root_cld)
plot_summary <- InF %>%
dplyr::select(trt, location, stage, row, shoot_wt, root_wt_inf) %>%
dplyr::group_by(trt, location, stage) %>%
dplyr:: summarise(
mean_root = mean(root_wt_inf, na.rm = TRUE),
se_root = sd(root_wt_inf, na.rm = TRUE) / sqrt(sum(!is.na(root_wt_inf))),
mean_shoot = mean(shoot_wt, na.rm = TRUE),
se_shoot = sd(shoot_wt, na.rm = TRUE) / sqrt(sum(!is.na(shoot_wt))),
.groups = "drop"
) %>%
dplyr:: left_join(shoot_cld_df, by = c("trt", "location", "stage")) %>%
dplyr:: left_join(root_cld_df, by = c("trt", "location", "stage"))
plotdata <- plot_summary %>%
dplyr:: select(trt, location, stage, mean_root, mean_shoot) %>%
pivot_longer(
cols = c(mean_root, mean_shoot),
names_to = "pool",
values_to = "amount"
) %>%
dplyr:: mutate(
pool = factor(pool,
levels = c("mean_root", "mean_shoot"),
labels = c("Root", "Shoot"))
)
# shoot error bars: centered on top of shoot segment
shoot_err <- plot_summary %>%
dplyr::mutate(
ymin = mean_shoot - se_shoot,
ymax = mean_shoot + se_shoot
)
# root error bars: root segment sits above shoot
root_err <- plot_summary %>%
dplyr::mutate(
root_top = mean_shoot + mean_root,
ymin = root_top - se_root,
ymax = root_top + se_root
)
label_df <- plot_summary %>%
dplyr::group_by(location, stage) %>%
dplyr::mutate(
panel_max = max(mean_shoot + mean_root + se_root, na.rm = TRUE),
shoot_lab_y = mean_shoot + se_shoot + 0.06 * panel_max,
root_lab_y = mean_shoot + mean_root + se_root + 0.06 * panel_max
) %>%
ungroup()
dodge_amt <- 0.09
p1 <- ggplot(plotdata, aes(x = factor(trt), y = amount, fill = pool)) +
geom_col(width = 0.60) +
geom_errorbar(
data = shoot_err,
aes(
x = as.numeric(factor(trt)) - dodge_amt,
ymin = ymin,
ymax = ymax
),
width = 0.12,
linewidth = 0.5,
inherit.aes = FALSE
) +
geom_errorbar(
data = root_err,
aes(
x = as.numeric(factor(trt)) + dodge_amt,
ymin = ymin,
ymax = ymax
),
width = 0.12,
linewidth = 0.5,
color= "grey40",
inherit.aes = FALSE
) +
geom_text(
data = label_df,
aes(
x = as.numeric(factor(trt)) - dodge_amt,
y = shoot_lab_y,
label = shoot_cld
),
inherit.aes = FALSE,
size = 5
) +
geom_text(
data = label_df,
aes(
x = as.numeric(factor(trt)) + dodge_amt,
y = root_lab_y,
label = toupper(root_cld)
),
inherit.aes = FALSE,
size = 3.8,
color= "grey40"
) +
facet_grid(vars(stage), vars(location), scales = "free_y") +
labs(
x = "Fertilizer-N treatment",
y = expression("Crop biomass (kg ha"^{-1}*")"),
fill = "Component"
) +
scale_fill_manual(
values = c(
"Root" = "#D8A24A",
"Shoot" = "#4A9B83"
)
)+
theme_Manjot() +
theme(axis.text.x = element_text(angle = 45, hjust = 1))
p1
#ggsave("cropbiomass.png", p1, device = ragg::agg_png, width = 12, height = 8, units = "in", dpi = 300)Shoot N uptake (Soil vs Fertilizer) graph
#Fertilizer-N in shoot models
model_v3_neb_shoot <- lm(FertN_shoot_kg_ha ~ trt + row, data = dat_Neb_V3)
model_v6_neb_shoot <- lm(FertN_shoot_kg_ha ~ trt + row, data = dat_Neb_V6)
model_v3_nh_shoot <- lm(FertN_shoot_kg_ha ~ trt + row, data = dat_NH_V3)
model_v6_nh_shoot <- lm(FertN_shoot_kg_ha ~ trt + row, data = dat_NH_V6)
emm_v3_neb_shoot <- emmeans(model_v3_neb_shoot, ~ trt)
emm_v6_neb_shoot <- emmeans(model_v6_neb_shoot, ~ trt)
emm_v3_nh_shoot <- emmeans(model_v3_nh_shoot, ~ trt)
emm_v6_nh_shoot <- emmeans(model_v6_nh_shoot, ~ trt)
cld_v3_neb_shoot <- cld(emm_v3_neb_shoot, Letters = letters, alpha = 0.05, sort= TRUE, reversed= TRUE)
cld_v6_neb_shoot <- cld(emm_v6_neb_shoot, Letters = letters, alpha = 0.05, sort= TRUE, reversed= TRUE)
cld_v3_nh_shoot <- cld(emm_v3_nh_shoot, Letters = letters, alpha = 0.05, sort= TRUE, reversed= TRUE)
cld_v6_nh_shoot <- cld(emm_v6_nh_shoot, Letters = letters, alpha = 0.05, sort= TRUE, reversed= TRUE)
shoot_cld_df <- bind_rows(
as.data.frame(cld_v3_neb_shoot) %>% mutate(stage = "V3", location = "NE"),
as.data.frame(cld_v6_neb_shoot) %>% mutate(stage = "V6", location = "NE"),
as.data.frame(cld_v3_nh_shoot) %>% mutate(stage = "V3", location = "NH"),
as.data.frame(cld_v6_nh_shoot) %>% mutate(stage = "V6", location = "NH")
) %>%
dplyr::rename(shoot_cld = .group) %>%
dplyr::mutate(shoot_cld = gsub(" ", "", shoot_cld)) %>%
dplyr::select(trt, location, stage, shoot_cld)
# SOM derived shoot N uptake models
model_v3_neb_som <- lm(Ndf_notfert ~ trt + row, data = dat_Neb_V3)
model_v6_neb_som <- lm(Ndf_notfert ~ trt + row, data = dat_Neb_V6)
model_v3_nh_som <- lm(Ndf_notfert ~ trt + row, data = dat_NH_V3)
model_v6_nh_som <- lm(Ndf_notfert ~ trt + row, data = dat_NH_V6)
#none of them significant by trt
anova(model_v3_neb_som)## Analysis of Variance Table
##
## Response: Ndf_notfert
## Df Sum Sq Mean Sq F value Pr(>F)
## trt 5 1.61753 0.32351 2.2357 0.1044
## row 3 0.21087 0.07029 0.4858 0.6972
## Residuals 15 2.17054 0.14470anova(model_v6_neb_som)## Analysis of Variance Table
##
## Response: Ndf_notfert
## Df Sum Sq Mean Sq F value Pr(>F)
## trt 5 124.188 24.838 1.5131 0.2443
## row 3 35.878 11.959 0.7285 0.5507
## Residuals 15 246.231 16.415anova(model_v3_nh_som)## Analysis of Variance Table
##
## Response: Ndf_notfert
## Df Sum Sq Mean Sq F value Pr(>F)
## trt 5 1.20499 0.24100 1.8566 0.1622
## row 3 0.59304 0.19768 1.5229 0.2493
## Residuals 15 1.94708 0.12981anova(model_v6_nh_som)## Analysis of Variance Table
##
## Response: Ndf_notfert
## Df Sum Sq Mean Sq F value Pr(>F)
## trt 5 76.897 15.379 1.5819 0.2288
## row 3 73.044 24.348 2.5044 0.1016
## Residuals 14 136.108 9.722emm_v3_neb_som <- emmeans(model_v3_neb_som, ~ trt)
emm_v6_neb_som <- emmeans(model_v6_neb_som, ~ trt)
emm_v3_nh_som <- emmeans(model_v3_nh_som, ~ trt)
emm_v6_nh_som <- emmeans(model_v6_nh_som, ~ trt)
cld_v3_neb_som <- cld(emm_v3_neb_som, Letters = letters, alpha = 0.05, sort= TRUE, reversed= TRUE)
cld_v6_neb_som <- cld(emm_v6_neb_som, Letters = letters, alpha = 0.05, sort= TRUE, reversed= TRUE)
cld_v3_nh_som <- cld(emm_v3_nh_som, Letters = letters, alpha = 0.05, sort= TRUE, reversed= TRUE)
cld_v6_nh_som <- cld(emm_v6_nh_som, Letters = letters, alpha = 0.05, sort= TRUE, reversed= TRUE)
som_cld_df <- bind_rows(
as.data.frame(cld_v3_neb_som) %>% mutate(stage = "V3", location = "NE"),
as.data.frame(cld_v6_neb_som) %>% mutate(stage = "V6", location = "NE"),
as.data.frame(cld_v3_nh_som) %>% mutate(stage = "V3", location = "NH"),
as.data.frame(cld_v6_nh_som) %>% mutate(stage = "V6", location = "NH")
) %>%
dplyr:: rename(som_cld = .group) %>%
dplyr:: mutate(som_cld = gsub(" ", "", som_cld)) %>%
dplyr:: select(trt, location, stage, som_cld)
plot_summary <- data_new %>%
dplyr::select(trt, location, stage, row, FertN_shoot_kg_ha, Ndf_notfert) %>%
dplyr::group_by(trt, location, stage) %>%
dplyr:: summarise(
mean_SOM = mean(Ndf_notfert, na.rm = TRUE),
se_SOM = sd(Ndf_notfert, na.rm = TRUE) / sqrt(sum(!is.na(Ndf_notfert))),
mean_shoot = mean(FertN_shoot_kg_ha, na.rm = TRUE),
se_shoot = sd(FertN_shoot_kg_ha, na.rm = TRUE) / sqrt(sum(!is.na(FertN_shoot_kg_ha))),
.groups = "drop"
) %>%
left_join(shoot_cld_df, by = c("trt", "location", "stage")) %>%
left_join(som_cld_df, by = c("trt", "location", "stage"))
plotdata <- plot_summary %>%
dplyr:: select(trt, location, stage, mean_SOM, mean_shoot) %>%
pivot_longer(
cols = c(mean_shoot, mean_SOM),
names_to = "pool",
values_to = "amount"
) %>%
dplyr:: mutate(
pool = factor(pool,
levels = c("mean_shoot", "mean_SOM"),
labels = c("Fertilizer-N", "Soil-derived N"))
)
# shoot error bars: centered on top of shoot segment
SOM_err <- plot_summary %>%
dplyr::mutate(
ymin = mean_SOM - se_SOM,
ymax = mean_SOM + se_SOM
)
# root error bars: root segment sits above shoot
shoot_err <- plot_summary %>%
dplyr::mutate(
shoot_top = mean_shoot + mean_SOM,
ymin = shoot_top - se_shoot,
ymax = shoot_top + se_shoot
)
label_df <- plot_summary %>%
dplyr::group_by(location, stage) %>%
dplyr::mutate(
panel_max = max(mean_shoot + mean_SOM + se_shoot, na.rm = TRUE),
SOM_lab_y = mean_SOM + se_SOM + 0.06 * panel_max,
shoot_lab_y = mean_shoot + mean_SOM + se_shoot + 0.06 * panel_max
) %>%
ungroup()
dodge_amt <- 0.08
p1<-ggplot(plotdata, aes(x = factor(trt), y = amount, fill = pool)) +
geom_col(width = 0.60) +
geom_errorbar(
data = SOM_err,
aes(
x = as.numeric(factor(trt)) + dodge_amt,
ymin = ymin,
ymax = ymax
),
width = 0.12,
linewidth = 0.5,
inherit.aes = FALSE
) +
geom_errorbar(
data = shoot_err,
aes(
x = as.numeric(factor(trt)) - dodge_amt,
ymin = ymin,
ymax = ymax
),
width = 0.12,
linewidth = 0.5,
inherit.aes = FALSE
) +
geom_text(
data = label_df,
aes(
x = as.numeric(factor(trt)) - dodge_amt,
y = shoot_lab_y,
label = shoot_cld
),
inherit.aes = FALSE,
size = 5
) +
facet_grid(vars(stage), vars(location), scales = "free_y") +
labs(
x = "Fertilizer-N treatment",
y = expression("Shoot N uptake (kg N ha"^{-1}*")"),
fill = "Source of N"
) +
scale_fill_manual(
values = c(
"Fertilizer-N" = "#68AEDF", # blue
"Soil-derived N" = "#CC9352" # medium gray
)
)+
theme_Manjot() +
theme(axis.text.x = element_text(angle = 45, hjust = 1))
p1
#ggsave("shootNuptakeexp2.png", p1, device = ragg::agg_png, width = 12, height = 8, units = "in", dpi = 300)Regression shoot N uptake graph
# Reshape data
plotdata_long <- dat_V6 %>%
pivot_longer(cols = c(Ndf_notfert, FertN_shoot_kg_ha),
names_to = "source",
values_to = "shoot_N_uptake") %>%
dplyr::mutate(
source = recode(source,
"Ndf_notfert" = "Soil-derived N",
"FertN_shoot_kg_ha" = "Fertilizer-N")
)
# Plot with per-facet regressions handled by ggpmisc
p1<-ggplot(plotdata_long, aes(x = totN_applied_kgha,
y = shoot_N_uptake,
color = source, shape=source)) +
# Linear for SOM-N
geom_smooth(
data = subset(plotdata_long, source == "Soil-derived N"),
aes(group = interaction(source, location)),
method = "lm", formula = y ~ x, se = TRUE, alpha=0.20, fill= "#CC9352"
) +
# Linear for Fertilizer-N
geom_smooth(
data = subset(plotdata_long, source == "Fertilizer-N"),
aes(group = interaction(source, location)),
method = "lm", formula = y ~ x, se = TRUE, alpha= 0.20, fill= "#68AEDF"
) +
geom_point(size=2.5, stroke=1) +
# Linear equations for SOM-N
stat_poly_eq(
data = subset(plotdata_long, source == "Soil-derived N"),
aes(group = interaction(source, location),
label = after_stat(
paste( rr.label, p.value.label, sep = "*', '*")
)),
formula = y ~ x, parse = TRUE,
geom = "label_npc", label.x = "left", label.y = 0.95,
vstep = 0.12, size = 4.5
) +
# Linear equations for Fertilizer-N
stat_poly_eq(
data = subset(plotdata_long, source == "Fertilizer-N"),
aes(group = interaction(source, location),
label = after_stat(
paste( rr.label, p.value.label, sep = "*', '*")
)),
formula = y ~ x, parse = TRUE,
geom = "label_npc", label.x = "left", label.y = 0.88,
vstep = 0.12, size = 4.5
) +
facet_wrap(~ location, ncol = 2,
labeller = as_labeller(c(NE = "NE", NH = "NH")))+
scale_color_manual(values = c("Soil-derived N" = "#CC9352",
"Fertilizer-N" = "#6BAEDF")) +
scale_fill_manual(
values = c(
"Soil-derived N" = "#CC9352",
"Fertilizer-N" = "#6BAEDF"
))+
scale_shape_manual(values = c("Soil-derived N" = 16,
"Fertilizer-N" = 17)) +
labs(
x = expression("Fertilizer Rate (kg N " ~ha^{-1} ~ ")"),
y = expression("Shoot N Uptake (kg N " ~ha^{-1} ~ ")"),
color = "Source of N", shape= "Source of N"
) +
theme_Manjot()+
theme(
axis.text = element_text(size = 15)
)
p1## Warning: Removed 49 rows containing non-finite outside the scale range
## (`stat_smooth()`).
## Removed 49 rows containing non-finite outside the scale range
## (`stat_smooth()`).## Warning: Removed 49 rows containing non-finite outside the scale range
## (`stat_poly_eq()`).
## Removed 49 rows containing non-finite outside the scale range
## (`stat_poly_eq()`).## Warning: No shared levels found between `names(values)` of the manual scale and the
## data's fill values.## Warning: Removed 98 rows containing missing values or values outside the scale range
## (`geom_point()`).
ggsave("NsourceAexp2.jpg", p1, device = ragg::agg_png, width = 12, height = 8, units = "in", dpi = 300)## Warning: Removed 49 rows containing non-finite outside the scale range
## (`stat_smooth()`).## Warning: Removed 49 rows containing non-finite outside the scale range
## (`stat_smooth()`).## Warning: Removed 49 rows containing non-finite outside the scale range
## (`stat_poly_eq()`).
## Removed 49 rows containing non-finite outside the scale range
## (`stat_poly_eq()`).## Warning: No shared levels found between `names(values)` of the manual scale and the
## data's fill values.## Warning: Removed 98 rows containing missing values or values outside the scale range
## (`geom_point()`).15N recovery graph
plotdata <- InF %>%
dplyr::select(
trt, location, stage, row,
per_Recov_wholesoil_tot,
perRecov_shoot,
perRecov_root
) %>%
dplyr:: group_by(trt, location, stage) %>%
dplyr:: summarise(
Root = mean(perRecov_root, na.rm = TRUE),
Shoot = mean(perRecov_shoot, na.rm = TRUE),
`Whole soil` = mean(per_Recov_wholesoil_tot, na.rm = TRUE),
.groups = "drop"
) %>%
pivot_longer(
cols = c(Root, Shoot, `Whole soil`),
names_to = "pool",
values_to = "amount"
) %>%
dplyr:: mutate(
pool = factor(
pool,
levels = c("Root", "Shoot", "Whole soil")
)
)
p1 <- ggplot(
plotdata,
aes(x = factor(trt), y = amount, fill = pool)
) +
geom_col(
position = "stack",
width = 0.60
) +
facet_grid(
vars(stage),
vars(location),
scales = "free_y"
) +
labs(
x = "Fertilizer treatment",
y = expression({}^{15}*"N recovery (%)"),
fill = "Pools"
) +
scale_fill_manual(
values = c(
"Root" = "#D8A24A",
"Shoot" = "#4A9B83",
"Whole soil" = "#8C78A8"
)
) +
theme_Manjot() +
theme(
axis.text = element_text(size = 15),
axis.text.x = element_text(
angle = 45,
hjust = 1
),
axis.title = element_text(size = 16)
)
p1
ggsave("15Nrecovery_npj.jpg", p1, device = ragg::agg_png, width = 10, height = 8, units = "in", dpi = 300)15N recovery in soil pools graph
plotdata <- data_new %>%
dplyr::select(
trt, location, stage, row,
per_Recov_MAOM_tot,
per_Recov_TDN_tot,
per_Recov_MBN_tot
) %>%
dplyr:: group_by(trt, location, stage) %>%
dplyr:: summarise(
mean_maomrec = mean(per_Recov_MAOM_tot, na.rm = TRUE),
mean_mbnrec = mean(per_Recov_MBN_tot, na.rm = TRUE),
mean_tdnrec = mean(per_Recov_TDN_tot, na.rm = TRUE),
.groups = "drop"
) %>%
pivot_longer(
cols = c(mean_maomrec, mean_mbnrec, mean_tdnrec),
names_to = "pool",
values_to = "amount"
) %>%
dplyr:: mutate(
pool = factor(
pool,
levels = c(
"mean_maomrec",
"mean_mbnrec",
"mean_tdnrec"
)
)
)
axis_caps <- plotdata %>%
dplyr::distinct(stage, location) %>%
dplyr::mutate(
trt = dplyr::first(plotdata$trt),
amount = dplyr::if_else(stage == "V3", 110, 60)
)
p1 <- ggplot(
plotdata,
aes(
x = factor(trt),
y = amount,
fill = pool
)
) +
geom_col(
position = "stack",
width = 0.60,
linewidth = 0.25
) +
facet_grid(
vars(stage),
vars(location),
scales = "free_y"
) +
scale_fill_manual(
values = c(
"mean_maomrec" = "#8E5A8A", # warm purple
"mean_mbnrec" = "#C48FB3", # light warm purple
"mean_tdnrec" = "#6678B8" # cool blue-purple
),
breaks = c(
"mean_maomrec",
"mean_mbnrec",
"mean_tdnrec"
),
labels = c(
"MAOM",
"MBN",
"TDN"
)
) +
# Invisible points force V3 to 90 and V6 to 60
geom_blank(
data = axis_caps,
aes(
x = factor(trt),
y = amount
),
inherit.aes = FALSE
) +
labs(
x = "Fertilizer treatment",
y = expression({}^{15}*"N recovery (%)"),
fill = "Soil pools")+
theme_Manjot() +
theme(
axis.text = element_text(size = 15),
axis.text.x = element_text(
angle = 45,
hjust = 1
),
axis.title = element_text(size = 16),
legend.title = element_text(size = 15),
legend.text = element_text(size = 14)
)
p1
#ggsave("recoverysoilpool2.jpg", p1, device = ragg::agg_png, width = 10, height = 8, units = "in", dpi = 300)Ndff in MAOM graph
#MAOM band model
model_v3_neb_maom <- lm(kgN_dff_MAOM_ha ~ trt + row, data = dat_Neb_V3_Band)
model_v6_neb_maom <- lm(kgN_dff_MAOM_ha~ trt + row, data = dat_Neb_V6_Band)
model_v3_nh_maom <- lm(kgN_dff_MAOM_ha~ trt + row, data = dat_NH_V3_Band)
model_v6_nh_maom <- lm(kgN_dff_MAOM_ha~ trt + row, data = dat_NH_V6_Band)
emm_v3_neb_maom <- emmeans(model_v3_neb_maom, ~ trt)
emm_v6_neb_maom <- emmeans(model_v6_neb_maom, ~ trt)
emm_v3_nh_maom <- emmeans(model_v3_nh_maom, ~ trt)
emm_v6_nh_maom <- emmeans(model_v6_nh_maom, ~ trt)
cld_v3_neb_maom <- cld(emm_v3_neb_maom, Letters = letters, alpha = 0.05, sort= TRUE, reversed= TRUE)
cld_v6_neb_maom <- cld(emm_v6_neb_maom, Letters = letters, alpha = 0.05, sort= TRUE, reversed= TRUE)
cld_v3_nh_maom <- cld(emm_v3_nh_maom, Letters = letters, alpha = 0.05, sort= TRUE, reversed= TRUE)
cld_v6_nh_maom <- cld(emm_v6_nh_maom, Letters = letters, alpha = 0.05, sort= TRUE, reversed= TRUE)
#-----------------------------
# MAOM furrow models
#-----------------------------
model_v3_neb_furrow <- lm(kgN_dff_MAOM_ha ~ trt + row, data= dat_Neb_V3_Inf)
model_v6_neb_furrow <- lm(kgN_dff_MAOM_ha~ trt + row, data = dat_Neb_V6_Inf)
model_v3_nh_furrow <- lm(kgN_dff_MAOM_ha ~ trt + row, data = dat_NH_V3_Inf)
model_v6_nh_furrow <- lm(kgN_dff_MAOM_ha ~ trt + row, data = dat_NH_V6_Inf)
emm_v3_neb_furrow <- emmeans(model_v3_neb_furrow, ~ trt)
emm_v6_neb_furrow <- emmeans(model_v6_neb_furrow, ~ trt)
emm_v3_nh_furrow <- emmeans(model_v3_nh_furrow, ~ trt)
emm_v6_nh_furrow <- emmeans(model_v6_nh_furrow, ~ trt)
cld_v3_neb_furrow <- cld(emm_v3_neb_furrow, Letters = letters, alpha = 0.05, sort= TRUE, reversed= TRUE)
cld_v6_neb_furrow <- cld(emm_v6_neb_furrow, Letters = letters, alpha = 0.05, sort= TRUE, reversed= TRUE)
cld_v3_nh_furrow <- cld(emm_v3_nh_furrow, Letters = letters, alpha = 0.05, sort= TRUE, reversed= TRUE)
cld_v6_nh_furrow <- cld(emm_v6_nh_furrow, Letters = letters, alpha = 0.05, sort= TRUE, reversed= TRUE)
# Combine significance letters for Band and Furrow
cld_df <- bind_rows(
as.data.frame(cld_v3_neb_maom) %>%
dplyr::mutate(stage = "V3", location = "NE", position = "Band"),
as.data.frame(cld_v6_neb_maom) %>%
dplyr::mutate(stage = "V6", location = "NE", position = "Band"),
as.data.frame(cld_v3_nh_maom) %>%
dplyr::mutate(stage = "V3", location = "NH", position = "Band"),
as.data.frame(cld_v6_nh_maom) %>%
dplyr::mutate(stage = "V6", location = "NH", position = "Band"),
as.data.frame(cld_v3_neb_furrow) %>%
dplyr::mutate(stage = "V3", location = "NE", position = "Furrow"),
as.data.frame(cld_v6_neb_furrow) %>%
dplyr:: mutate(stage = "V6", location = "NE", position = "Furrow"),
as.data.frame(cld_v3_nh_furrow) %>%
dplyr::mutate(stage = "V3", location = "NH", position = "Furrow"),
as.data.frame(cld_v6_nh_furrow) %>%
dplyr:: mutate(stage = "V6", location = "NH", position = "Furrow")
) %>%
dplyr::transmute(
trt = as.character(trt),
location,
stage,
position,
cld = gsub("\\s+", "", .group)
)
# Treatment order
trt_levels <- c(
"InF-5", "InF-10", "UAN-74",
"Both-79", "UAN-148", "Both-158"
)
# Standardize position labels in CLD data
cld_df <- cld_df %>%
dplyr::mutate(
trt = trimws(as.character(trt)),
location = as.character(location),
stage = as.character(stage),
position = trimws(as.character(position))
)
# Summarize MAOM by treatment, location, stage, and position
plotdata <- data_new %>%
dplyr::mutate(
trt = trimws(as.character(trt)),
location = as.character(location),
stage = as.character(stage),
position = dplyr::case_when(
trimws(as.character(position)) %in%
c("In-row", "Furrow", "InF") ~ "Furrow",
trimws(as.character(position)) == "Band" ~ "Band",
TRUE ~ NA_character_
)
) %>%
dplyr::filter(
!is.na(position),
!is.na(kgN_dff_MAOM_ha)
) %>%
dplyr::group_by(trt, location, stage, position) %>%
dplyr::summarise(
mean_maom = mean(kgN_dff_MAOM_ha),
se_maom = if (dplyr::n() > 1) {
sd(kgN_dff_MAOM_ha) / sqrt(dplyr::n())
} else {
NA_real_
},
.groups = "drop"
) %>%
# Join significance letters AFTER standardizing position
dplyr::left_join(
cld_df,
by = c("trt", "location", "stage", "position")
) %>%
dplyr::mutate(
trt = factor(trt, levels = trt_levels),
position = factor(
position,
levels = c("Band", "Furrow")
)
) %>%
# Calculate significance letter heights
dplyr::group_by(location, stage) %>%
dplyr::mutate(
panel_max = max(
mean_maom + dplyr::coalesce(se_maom, 0),
na.rm = TRUE
),
label_y = mean_maom +
dplyr::coalesce(se_maom, 0) +
dplyr::if_else(
position == "Furrow",
0.10 * max(panel_max, 1),
0.06 * max(panel_max, 1)
)
) %>%
dplyr::ungroup()
# Error bars
error_df <- plotdata %>%
dplyr::mutate(
ymin = mean_maom - se_maom,
ymax = mean_maom + se_maom
)
# Significance letters
label_long <- plotdata %>%
dplyr::select(
trt, location, stage, position,
label_y, cld
)
# Identical dodge for all layers
dodge_position <- ggplot2::position_dodge(width = 0.75)
# Plot
p1 <- ggplot2::ggplot(
plotdata,
ggplot2::aes(
x = trt,
y = mean_maom,
fill = position,
group = position
)
) +
# Bars
ggplot2::geom_col(
position = dodge_position,
width = 0.60
) +
# Error bars
ggplot2::geom_errorbar(
data = error_df,
ggplot2::aes(
x = trt,
ymin = ymin,
ymax = ymax,
group = position
),
position = dodge_position,
width = 0.08,
linewidth = 0.5,
color = "grey20",
inherit.aes = FALSE,
na.rm = TRUE
) +
# Significance letters
ggplot2::geom_text(
data = label_long,
ggplot2::aes(
x = trt,
y = label_y,
label = cld,
group = position
),
position = dodge_position,
inherit.aes = FALSE,
size = 5,
color = "grey15",
na.rm = TRUE
) +
ggplot2::facet_grid(
ggplot2::vars(stage),
ggplot2::vars(location),
scales = "free_y"
) +
ggplot2::scale_fill_manual(
values = c(
"Band" = "#8E5A8A",
"Furrow" = "#6678B8"
),
breaks = c("Band", "Furrow")
) +
ggplot2::labs(
x = "Fertilizer-N treatment",
y = expression(
"Amount Ndff in MAOM (kg N ha"^{-1} * ")"
),
fill = "Position"
) +
theme_Manjot() +
ggplot2::theme(
axis.text = ggplot2::element_text(size = 15),
axis.text.x = ggplot2::element_text(
angle = 45,
hjust = 1
),
axis.title = ggplot2::element_text(size = 16),
legend.title = ggplot2::element_text(size = 16),
legend.text = ggplot2::element_text(size = 15)
)
p1
ggsave("MAOMdffexp2.jpg", p1, device = ragg::agg_png, width = 12, height = 8, units = "in", dpi = 300)% recovery in MAOM graph
# Fit separate models for each site x stage
#-----------------------------
model_v3_neb_maom <- lm(per_Recov_MAOM_tot ~ trt + row, data = dat_Neb_V3)
model_v6_neb_maom <- lm(per_Recov_MAOM_tot ~ trt + row, data = dat_Neb_V6)
model_v3_nh_maom <- lm(per_Recov_MAOM_tot ~ trt + row, data = dat_NH_V3)
model_v6_nh_maom <- lm(per_Recov_MAOM_tot ~ trt + row, data = dat_NH_V6)
#-----------------------------
# emmeans: compare trt
#-----------------------------
emm_v3_neb_maom <- emmeans(model_v3_neb_maom, ~ trt )
emm_v6_neb_maom <- emmeans(model_v6_neb_maom, ~ trt )
emm_v3_nh_maom <- emmeans(model_v3_nh_maom, ~ trt )
emm_v6_nh_maom <- emmeans(model_v6_nh_maom, ~ trt )
#-----------------------------
# CLD letters
#-----------------------------
cld_v3_neb_maom <- cld(emm_v3_neb_maom, Letters = letters, alpha = 0.05, sort=TRUE, reversed=TRUE)
cld_v6_neb_maom <- cld(emm_v6_neb_maom, Letters = letters, alpha = 0.05,sort=TRUE, reversed=TRUE)
cld_v3_nh_maom <- cld(emm_v3_nh_maom, Letters = letters, alpha = 0.05,sort=TRUE, reversed=TRUE)
cld_v6_nh_maom <- cld(emm_v6_nh_maom, Letters = letters, alpha = 0.05,sort=TRUE, reversed=TRUE)
#-----------------------------
# Convert CLD results to dataframes
#-----------------------------
maom_cld_df <- bind_rows(
as.data.frame(cld_v3_neb_maom) %>% dplyr::mutate(stage = "V3", location = "NE"),
as.data.frame(cld_v6_neb_maom) %>% dplyr::mutate(stage = "V6", location = "NE"),
as.data.frame(cld_v3_nh_maom) %>% dplyr::mutate(stage = "V3", location = "NH"),
as.data.frame(cld_v6_nh_maom) %>% dplyr::mutate(stage = "V6", location = "NH")
) %>%
dplyr::rename(maom_cld = .group) %>%
dplyr::mutate(maom_cld = gsub(" ", "", maom_cld)) %>%
dplyr::select(trt, location, stage, maom_cld)
label_df <- data_new %>%
dplyr::filter(stage %in% c("V3", "V6")) %>%
dplyr::group_by(location, stage, trt) %>%
dplyr::summarise(
y_pos = boxplot.stats(per_Recov_MAOM_tot)$stats[5], # upper whisker
.groups = "drop"
) %>%
left_join(maom_cld_df, by = c("trt", "location", "stage")) %>%
dplyr::group_by(location, stage) %>%
dplyr::mutate(
panel_max = max(y_pos, na.rm = TRUE),
y_pos = y_pos + 0.1 * panel_max
) %>%
ungroup()
ggplot(data_new %>% dplyr::filter(stage %in% c("V3", "V6")),
aes(x = trt, y = per_Recov_MAOM_tot)) +
geom_boxplot(width = 0.6) +
facet_grid(vars(stage), vars(location), scales = "free") +
geom_text(
data = label_df,
aes(x = trt, y = y_pos, label = maom_cld),
position = position_dodge(width = 0.8),
inherit.aes = FALSE,
size = 5
) +
theme_Manjot() +
labs(
x = "Fertilizer treatment",
y = "% recovery in MAOM",
title = "Recovery in MAOM"
) +
theme(axis.text.x = element_text(angle = 45, hjust = 1))## Warning: Removed 96 rows containing non-finite outside the scale range
## (`stat_boxplot()`).
## Total N in crop biomass graph
model_v3_neb_shoot <- lm(totN_shoot_kgha ~ trt + row, data = dat_Neb_V3)
model_v6_neb_shoot <- lm(totN_shoot_kgha ~ trt + row, data = dat_Neb_V6)
model_v3_nh_shoot <- lm(totN_shoot_kgha ~ trt + row, data = dat_NH_V3)
model_v6_nh_shoot <- lm(totN_shoot_kgha ~ trt + row, data = dat_NH_V6)
emm_v3_neb_shoot <- emmeans(model_v3_neb_shoot, ~ trt)
emm_v6_neb_shoot <- emmeans(model_v6_neb_shoot, ~ trt)
emm_v3_nh_shoot <- emmeans(model_v3_nh_shoot, ~ trt)
emm_v6_nh_shoot <- emmeans(model_v6_nh_shoot, ~ trt)
cld_v3_neb_shoot <- cld(emm_v3_neb_shoot, Letters = letters, alpha = 0.05, sort= TRUE, reversed= TRUE)
cld_v6_neb_shoot <- cld(emm_v6_neb_shoot, Letters = letters, alpha = 0.05, sort= TRUE, reversed= TRUE)
cld_v3_nh_shoot <- cld(emm_v3_nh_shoot, Letters = letters, alpha = 0.05, sort= TRUE, reversed= TRUE)
cld_v6_nh_shoot <- cld(emm_v6_nh_shoot, Letters = letters, alpha = 0.05, sort= TRUE, reversed= TRUE)
shoot_cld_df <- bind_rows(
as.data.frame(cld_v3_neb_shoot) %>% mutate(stage = "V3", location = "NE"),
as.data.frame(cld_v6_neb_shoot) %>% mutate(stage = "V6", location = "NE"),
as.data.frame(cld_v3_nh_shoot) %>% mutate(stage = "V3", location = "NH"),
as.data.frame(cld_v6_nh_shoot) %>% mutate(stage = "V6", location = "NH")
) %>%
dplyr::rename(shoot_cld = .group) %>%
dplyr::mutate(shoot_cld = gsub(" ", "", shoot_cld)) %>%
dplyr::select(trt, location, stage, shoot_cld)
#-----------------------------
# ROOT models
#-----------------------------
model_v3_neb_root <- lm(totN_root_kgha ~ trt + row, data = dat_Neb_V3_Inf)
model_v6_neb_root <- lm(totN_root_kgha ~ trt + row, data = dat_Neb_V6_Inf)
model_v3_nh_root <- lm(totN_root_kgha ~ trt + row, data = dat_NH_V3_Inf)
model_v6_nh_root <- lm(totN_root_kgha ~ trt + row, data = dat_NH_V6_Inf)
emm_v3_neb_root <- emmeans(model_v3_neb_root, ~ trt)
emm_v6_neb_root <- emmeans(model_v6_neb_root, ~ trt)
emm_v3_nh_root <- emmeans(model_v3_nh_root, ~ trt)
emm_v6_nh_root <- emmeans(model_v6_nh_root, ~ trt)
cld_v3_neb_root <- cld(emm_v3_neb_root, Letters = letters, alpha = 0.05, sort= TRUE, reversed= TRUE)
cld_v6_neb_root <- cld(emm_v6_neb_root, Letters = letters, alpha = 0.05, sort= TRUE, reversed= TRUE)
cld_v3_nh_root <- cld(emm_v3_nh_root, Letters = letters, alpha = 0.05, sort= TRUE, reversed= TRUE)
cld_v6_nh_root <- cld(emm_v6_nh_root, Letters = letters, alpha = 0.05, sort= TRUE, reversed= TRUE)
root_cld_df <- bind_rows(
as.data.frame(cld_v3_neb_root) %>% mutate(stage = "V3", location = "NE"),
as.data.frame(cld_v6_neb_root) %>% mutate(stage = "V6", location = "NE"),
as.data.frame(cld_v3_nh_root) %>% mutate(stage = "V3", location = "NH"),
as.data.frame(cld_v6_nh_root) %>% mutate(stage = "V6", location = "NH")
) %>%
dplyr:: rename(root_cld = .group) %>%
dplyr:: mutate(root_cld = gsub(" ", "", root_cld)) %>%
dplyr:: select(trt, location, stage, root_cld)
plot_summary <- InF %>%
dplyr::select(trt, location, stage, row, totN_shoot_kgha, totN_root_kgha ) %>%
dplyr::group_by(trt, location, stage) %>%
dplyr:: summarise(
mean_root = mean(totN_root_kgha, na.rm = TRUE),
se_root = sd(totN_root_kgha, na.rm = TRUE) / sqrt(sum(!is.na(totN_root_kgha))),
mean_shoot = mean(totN_shoot_kgha, na.rm = TRUE),
se_shoot = sd(totN_shoot_kgha, na.rm = TRUE) / sqrt(sum(!is.na(totN_shoot_kgha))),
.groups = "drop"
) %>%
left_join(shoot_cld_df, by = c("trt", "location", "stage")) %>%
left_join(root_cld_df, by = c("trt", "location", "stage"))
plotdata <- plot_summary %>%
dplyr:: select(trt, location, stage, mean_root, mean_shoot) %>%
pivot_longer(
cols = c(mean_root, mean_shoot),
names_to = "pool",
values_to = "amount"
) %>%
dplyr:: mutate(
pool = factor(pool,
levels = c("mean_root", "mean_shoot"),
labels = c("Root", "Shoot"))
)
# shoot error bars: centered on top of shoot segment
shoot_err <- plot_summary %>%
dplyr::mutate(
ymin = mean_shoot - se_shoot,
ymax = mean_shoot + se_shoot
)
# root error bars: root segment sits above shoot
root_err <- plot_summary %>%
dplyr::mutate(
root_top = mean_shoot + mean_root,
ymin = root_top - se_root,
ymax = root_top + se_root
)
label_df <- plot_summary %>%
dplyr::group_by(location, stage) %>%
dplyr::mutate(
panel_max = max(mean_shoot + mean_root + se_root, na.rm = TRUE),
shoot_lab_y = mean_shoot + se_shoot + 0.06 * panel_max,
root_lab_y = mean_shoot + mean_root + se_root + 0.06 * panel_max
) %>%
ungroup()
p1 <- ggplot(plotdata, aes(x = factor(trt), y = amount, fill = pool)) +
geom_col(width = 0.60) +
geom_errorbar(
data = shoot_err,
aes(
x = as.numeric(factor(trt)) - dodge_amt,
ymin = ymin,
ymax = ymax
),
width = 0.12,
linewidth = 0.5,
inherit.aes = FALSE
) +
geom_errorbar(
data = root_err,
aes(
x = as.numeric(factor(trt)) + dodge_amt,
ymin = ymin,
ymax = ymax
),
width = 0.12,
linewidth = 0.5,
color= "grey40",
inherit.aes = FALSE
) +
geom_text(
data = label_df,
aes(
x = as.numeric(factor(trt)) - dodge_amt,
y = shoot_lab_y,
label = shoot_cld
),
inherit.aes = FALSE,
size = 5
) +
geom_text(
data = label_df,
aes(
x = as.numeric(factor(trt)) + dodge_amt,
y = root_lab_y,
label = toupper(root_cld)
),
inherit.aes = FALSE,
size = 3.8,
color= "grey40"
) +
facet_grid(vars(stage), vars(location), scales = "free_y") +
labs(
x = "Fertilizer-N treatment",
y = expression("Total N in crop biomass (kg ha"^{-1}*")"),
fill = "Component"
) +
scale_fill_manual(
values = c(
"Root" = "#D8A24A",
"Shoot" = "#4A9B83"
)
)+
theme_Manjot() +
theme(axis.text.x = element_text(angle = 45, hjust = 1))
p1
ggsave("TotNcropbiomass.png", p1, device = ragg::agg_png, width = 12, height = 8, units = "in", dpi = 300)Mean %15N recovery of all pools (table)
table_recovery <- InF%>%
pivot_longer(
cols = c(
perRecov_shoot,
perRecov_root,
per_Recov_MBN_tot,
per_Recov_MAOM_tot,
per_Recov_TDN_tot,
per_Recov_wholesoil_tot
),
names_to = "pool",
values_to = "recovery"
) %>%
dplyr::mutate(
pool = recode(pool,
perRecov_shoot = "Shoot",
perRecov_root = "Root",
per_Recov_MBN_tot = "MBN",
per_Recov_MAOM_tot= "MAOM",
per_Recov_TDN_tot = "TDN",
per_Recov_wholesoil_tot = "Whole soil"
)
) %>%
dplyr::group_by(trt, stage, location, pool) %>%
dplyr::summarise(
n = sum(!is.na(recovery)),
mean = mean(recovery, na.rm = TRUE),
se = sd(recovery, na.rm = TRUE) / sqrt(n),
.groups = "drop"
) %>%
dplyr:: mutate(
value = sprintf("%.1f (%.1f)", mean, se)
) %>%
dplyr::select(trt, stage, location, pool, value) %>%
pivot_wider(
names_from = pool,
values_from = value
)
View(table_recovery)
recovery_long <- InF %>%
tidyr::pivot_longer(
cols = c(
perRecov_shoot,
perRecov_root,
per_Recov_MBN_tot,
per_Recov_MAOM_tot,
per_Recov_TDN_tot,
per_Recov_wholesoil_tot
),
names_to = "pool",
values_to = "recovery"
) %>%
dplyr::mutate(
pool = dplyr::recode(
pool,
perRecov_shoot = "Shoot",
perRecov_root = "Root",
per_Recov_MBN_tot = "MBN",
per_Recov_MAOM_tot = "MAOM",
per_Recov_TDN_tot = "TDN",
per_Recov_wholesoil_tot = "Whole soil"
),
trt = factor(trt),
row = factor(row)
)
recovery_cld <- recovery_long %>%
dplyr::filter(
!is.na(recovery),
!is.na(trt),
!is.na(row)
) %>%
dplyr::group_by(location, stage, pool) %>%
dplyr::group_modify(~ {
model <- lm(
recovery ~ trt + row,
data = .x
)
emm <- emmeans::emmeans(model, ~ trt)
cld_result <- multcomp::cld(
emm,
Letters = letters,
alpha = 0.05,
sort = TRUE,
reversed = TRUE
)
data.frame(
trt = as.character(cld_result$trt),
cld = gsub("\\s+", "", cld_result$.group)
)
}) %>%
dplyr::ungroup()
table_recovery <- recovery_long %>%
dplyr::group_by(trt, stage, location, pool) %>%
dplyr::summarise(
n = sum(!is.na(recovery)),
mean = mean(recovery, na.rm = TRUE),
se = sd(recovery, na.rm = TRUE) / sqrt(n),
.groups = "drop"
) %>%
dplyr::mutate(
trt = as.character(trt)
) %>%
dplyr::left_join(
recovery_cld,
by = c("trt", "stage", "location", "pool")
) %>%
dplyr::mutate(
value = dplyr::case_when(
n > 1 ~ sprintf("%.1f (%.1f)", mean, se),
n == 1 ~ sprintf("%.1f (NA)", mean),
TRUE ~ NA_character_
),
value = dplyr::if_else(
!is.na(cld) & !is.na(value),
paste0(value, " ", cld),
value
)
) %>%
dplyr::select(
trt, stage, location, pool, value
) %>%
tidyr::pivot_wider(
names_from = pool,
values_from = value
)
table_recovery## # A tibble: 24 × 9
## trt stage location MAOM MBN Root Shoot TDN `Whole soil`
## <chr> <fct> <fct> <chr> <chr> <chr> <chr> <chr> <chr>
## 1 InF-5 V3 NE 23.4 (3.6) a 0.3 (0.1)… 3.1 … 13.3… 0.4 … 44.5 (20.5)…
## 2 InF-5 V3 NH 59.1 (11.5) a 0.4 (0.1)… 3.8 … 17.5… 0.3 … 87.5 (3.5) a
## 3 InF-5 V6 NE 11.3 (1.0) a 0.0 (0.0)… 1.5 … 22.0… 0.0 … 13.3 (1.6) a
## 4 InF-5 V6 NH 23.3 (2.8) a 0.1 (0.0)… 5.8 … 20.8… 0.1 … 35.3 (2.1) a
## 5 InF-10 V3 NE 15.7 (1.2) ab 0.1 (0.0)… 3.9 … 15.7… 0.5 … 42.0 (8.3) a
## 6 InF-10 V3 NH 27.8 (6.2) b 0.4 (0.2)… 3.6 … 14.9… 0.7 … 61.2 (2.7) b
## 7 InF-10 V6 NE 5.9 (0.5) ab 0.0 (0.0)… 2.5 … 24.0… 0.0 … 8.3 (1.0) a
## 8 InF-10 V6 NH 16.8 (6.4) ab 0.1 (0.0)… 5.7 … 24.9… 0.1 … 21.4 (3.2) a
## 9 UAN-74 V3 NE 9.5 (3.2) b 0.2 (0.2)… 0.2 … 1.0 … 7.3 … 37.3 (6.4) a
## 10 UAN-74 V3 NH 12.3 (1.7) b 3.5 (1.4)… 0.1 … 0.3 … 15.5… 56.2 (7.1) b
## # ℹ 14 more rowsknitr::kable(
table_recovery,
format = "html",
caption = "Fertilizer-N recovery (%) across soil and plant pools"
)| trt | stage | location | MAOM | MBN | Root | Shoot | TDN | Whole soil |
|---|---|---|---|---|---|---|---|---|
| InF-5 | V3 | NE | 23.4 (3.6) a | 0.3 (0.1) a | 3.1 (0.5) a | 13.3 (1.9) a | 0.4 (0.1) b | 44.5 (20.5) a |
| InF-5 | V3 | NH | 59.1 (11.5) a | 0.4 (0.1) a | 3.8 (0.7) a | 17.5 (3.5) a | 0.3 (0.1) d | 87.5 (3.5) a |
| InF-5 | V6 | NE | 11.3 (1.0) a | 0.0 (0.0) a | 1.5 (0.2) ab | 22.0 (6.4) a | 0.0 (0.0) b | 13.3 (1.6) a |
| InF-5 | V6 | NH | 23.3 (2.8) a | 0.1 (0.0) a | 5.8 (1.3) a | 20.8 (5.1) a | 0.1 (0.0) c | 35.3 (2.1) a |
| InF-10 | V3 | NE | 15.7 (1.2) ab | 0.1 (0.0) a | 3.9 (1.0) a | 15.7 (1.6) a | 0.5 (0.2) b | 42.0 (8.3) a |
| InF-10 | V3 | NH | 27.8 (6.2) b | 0.4 (0.2) a | 3.6 (0.3) a | 14.9 (1.3) a | 0.7 (0.2) d | 61.2 (2.7) b |
| InF-10 | V6 | NE | 5.9 (0.5) ab | 0.0 (0.0) a | 2.5 (0.4) a | 24.0 (3.5) a | 0.0 (0.0) b | 8.3 (1.0) a |
| InF-10 | V6 | NH | 16.8 (6.4) ab | 0.1 (0.0) a | 5.7 (0.8) a | 24.9 (3.6) a | 0.1 (0.0) c | 21.4 (3.2) a |
| UAN-74 | V3 | NE | 9.5 (3.2) b | 0.2 (0.2) a | 0.2 (0.1) b | 1.0 (0.5) b | 7.3 (2.9) ab | 37.3 (6.4) a |
| UAN-74 | V3 | NH | 12.3 (1.7) b | 3.5 (1.4) a | 0.1 (0.0) b | 0.3 (0.1) b | 15.5 (2.3) bc | 56.2 (7.1) b |
| UAN-74 | V6 | NE | 3.7 (0.5) b | 0.0 (0.0) a | 0.6 (0.2) bc | 23.4 (3.8) a | 0.1 (0.0) ab | 9.9 (1.0) a |
| UAN-74 | V6 | NH | 7.0 (1.7) b | 3.2 (3.2) a | 0.8 (0.3) b | 5.5 (2.9) b | 2.5 (0.8) bc | 25.2 (4.9) a |
| Both-79 | V3 | NE | 13.5 (2.3) ab | 1.5 (1.0) a | 0.4 (0.1) b | 1.8 (0.4) b | 12.4 (2.8) a | 30.4 (2.6) a |
| Both-79 | V3 | NH | 10.4 (1.6) b | 3.6 (2.1) a | 0.3 (0.1) b | 1.5 (0.4) b | 29.1 (5.9) a | 43.6 (4.7) b |
| Both-79 | V6 | NE | 5.2 (1.5) ab | 0.8 (0.6) a | 0.7 (0.1) bc | 15.1 (3.7) a | 2.6 (0.9) a | 14.1 (3.9) a |
| Both-79 | V6 | NH | 10.9 (1.2) ab | 0.6 (0.6) a | 0.6 (0.1) b | 5.1 (1.5) b | 17.5 (2.6) a | 35.6 (5.1) a |
| UAN-148 | V3 | NE | 9.8 (3.2) b | 0.3 (0.2) a | 0.2 (0.0) b | 0.7 (0.2) b | 4.3 (0.8) ab | 34.3 (7.2) a |
| UAN-148 | V3 | NH | 10.9 (1.3) b | 1.4 (0.9) a | 0.0 (0.0) b | 0.3 (0.1) b | 10.9 (3.0) cd | 47.1 (2.3) b |
| UAN-148 | V6 | NE | 5.2 (1.0) ab | 0.0 (0.0) a | 0.5 (0.1) c | 12.1 (1.3) a | 1.2 (1.1) ab | 15.3 (3.2) a |
| UAN-148 | V6 | NH | 10.6 (1.1) ab | 1.0 (0.7) a | 0.9 (0.2) b | 3.2 (0.9) b | 1.3 (0.4) bc | 35.3 (2.6) a |
| Both-158 | V3 | NE | 12.3 (2.0) ab | 0.8 (0.5) a | 0.4 (0.0) b | 1.5 (0.2) b | 13.6 (4.3) a | 34.2 (3.0) a |
| Both-158 | V3 | NH | 12.4 (1.8) b | 1.5 (1.0) a | 0.3 (0.0) b | 1.0 (0.2) b | 25.2 (2.9) ab | 36.7 (11.3) b |
| Both-158 | V6 | NE | 7.3 (2.5) ab | 0.0 (0.0) a | 0.7 (0.1) bc | 12.8 (1.6) a | 1.9 (0.5) ab | 16.0 (4.7) a |
| Both-158 | V6 | NH | 10.3 (1.0) ab | 0.4 (0.2) a | 1.2 (0.1) b | 6.4 (0.9) b | 5.1 (0.5) b | 28.0 (1.6) a |
N derived from fertilizer in crop biomass graph
model_v3_neb_shoot <- lm(FertN_shoot_kg_ha ~ trt + row, data = dat_Neb_V3)
model_v6_neb_shoot <- lm(FertN_shoot_kg_ha ~ trt + row, data = dat_Neb_V6)
model_v3_nh_shoot <- lm(FertN_shoot_kg_ha ~ trt + row, data = dat_NH_V3)
model_v6_nh_shoot <- lm(FertN_shoot_kg_ha ~ trt + row, data = dat_NH_V6)
emm_v3_neb_shoot <- emmeans(model_v3_neb_shoot, ~ trt)
emm_v6_neb_shoot <- emmeans(model_v6_neb_shoot, ~ trt)
emm_v3_nh_shoot <- emmeans(model_v3_nh_shoot, ~ trt)
emm_v6_nh_shoot <- emmeans(model_v6_nh_shoot, ~ trt)
cld_v3_neb_shoot <- cld(emm_v3_neb_shoot, Letters = letters, alpha = 0.05,sort= TRUE, reversed= TRUE)
cld_v6_neb_shoot <- cld(emm_v6_neb_shoot, Letters = letters, alpha = 0.05,sort= TRUE, reversed= TRUE)
cld_v3_nh_shoot <- cld(emm_v3_nh_shoot, Letters = letters, alpha = 0.05,sort= TRUE, reversed= TRUE)
cld_v6_nh_shoot <- cld(emm_v6_nh_shoot, Letters = letters, alpha = 0.05,sort= TRUE, reversed= TRUE)
shoot_cld_df <- bind_rows(
as.data.frame(cld_v3_neb_shoot) %>% mutate(stage = "V3", location = "NE"),
as.data.frame(cld_v6_neb_shoot) %>% mutate(stage = "V6", location = "NE"),
as.data.frame(cld_v3_nh_shoot) %>% mutate(stage = "V3", location = "NH"),
as.data.frame(cld_v6_nh_shoot) %>% mutate(stage = "V6", location = "NH")
) %>%
dplyr::rename(shoot_cld = .group) %>%
dplyr::mutate(shoot_cld = gsub(" ", "", shoot_cld)) %>%
dplyr::select(trt, location, stage, shoot_cld)
#-----------------------------
# ROOT models
#-----------------------------
model_v3_neb_root <- lm(FertN_root_kg_ha ~ trt + row, data = dat_Neb_V3_Inf)
model_v6_neb_root <- lm(FertN_root_kg_ha~ trt + row, data = dat_Neb_V6_Inf)
model_v3_nh_root <- lm(FertN_root_kg_ha ~ trt + row, data = dat_NH_V3_Inf)
model_v6_nh_root <- lm(FertN_root_kg_ha ~ trt + row, data = dat_NH_V6_Inf)
emm_v3_neb_root <- emmeans(model_v3_neb_root, ~ trt)
emm_v6_neb_root <- emmeans(model_v6_neb_root, ~ trt)
emm_v3_nh_root <- emmeans(model_v3_nh_root, ~ trt)
emm_v6_nh_root <- emmeans(model_v6_nh_root, ~ trt)
cld_v3_neb_root <- cld(emm_v3_neb_root, Letters = letters, alpha = 0.05,sort= TRUE, reversed= TRUE)
cld_v6_neb_root <- cld(emm_v6_neb_root, Letters = letters, alpha = 0.05,sort= TRUE, reversed= TRUE)
cld_v3_nh_root <- cld(emm_v3_nh_root, Letters = letters, alpha = 0.05,sort= TRUE, reversed= TRUE)
cld_v6_nh_root <- cld(emm_v6_nh_root, Letters = letters, alpha = 0.05,sort= TRUE, reversed= TRUE)
root_cld_df <- bind_rows(
as.data.frame(cld_v3_neb_root) %>% mutate(stage = "V3", location = "NE"),
as.data.frame(cld_v6_neb_root) %>% mutate(stage = "V6", location = "NE"),
as.data.frame(cld_v3_nh_root) %>% mutate(stage = "V3", location = "NH"),
as.data.frame(cld_v6_nh_root) %>% mutate(stage = "V6", location = "NH")
) %>%
dplyr:: rename(root_cld = .group) %>%
dplyr:: mutate(root_cld = gsub(" ", "", root_cld)) %>%
dplyr:: select(trt, location, stage, root_cld)
plot_summary <- InF %>%
dplyr::select(trt, location, stage, row, FertN_shoot_kg_ha, FertN_root_kg_ha) %>%
dplyr::group_by(trt, location, stage) %>%
dplyr:: summarise(
mean_root = mean(FertN_root_kg_ha, na.rm = TRUE),
se_root = sd(FertN_root_kg_ha, na.rm = TRUE) / sqrt(sum(!is.na(FertN_root_kg_ha))),
mean_shoot = mean(FertN_shoot_kg_ha, na.rm = TRUE),
se_shoot = sd(FertN_shoot_kg_ha, na.rm = TRUE) / sqrt(sum(!is.na(FertN_shoot_kg_ha))),
.groups = "drop"
) %>%
left_join(shoot_cld_df, by = c("trt", "location", "stage")) %>%
left_join(root_cld_df, by = c("trt", "location", "stage"))
plotdata <- plot_summary %>%
dplyr:: select(trt, location, stage, mean_root, mean_shoot) %>%
pivot_longer(
cols = c(mean_root, mean_shoot),
names_to = "pool",
values_to = "amount"
) %>%
dplyr:: mutate(
pool = factor(pool,
levels = c("mean_root", "mean_shoot"),
labels = c("Root", "Shoot"))
)
# shoot error bars: centered on top of shoot segment
shoot_err <- plot_summary %>%
dplyr::mutate(
ymin = mean_shoot - se_shoot,
ymax = mean_shoot + se_shoot
)
# root error bars: root segment sits above shoot
root_err <- plot_summary %>%
dplyr::mutate(
root_top = mean_shoot + mean_root,
ymin = root_top - se_root,
ymax = root_top + se_root
)
label_df <- plot_summary %>%
dplyr::group_by(location, stage) %>%
dplyr::mutate(
panel_max = max(mean_shoot + mean_root + se_root, na.rm = TRUE),
shoot_lab_y = mean_shoot + se_shoot + 0.06 * panel_max,
root_lab_y = mean_shoot + mean_root + se_root + 0.06 * panel_max
) %>%
ungroup()
p1 <- ggplot(plotdata, aes(x = factor(trt), y = amount, fill = pool)) +
geom_col(width = 0.60) +
geom_errorbar(
data = shoot_err,
aes(
x = as.numeric(factor(trt)) - 1.3*dodge_amt,
ymin = ymin,
ymax = ymax
),
width = 0.12,
linewidth = 0.5,
inherit.aes = FALSE
) +
geom_errorbar(
data = root_err,
aes(
x = as.numeric(factor(trt)) + 1.3*dodge_amt,
ymin = ymin,
ymax = ymax
),
width = 0.12,
linewidth = 0.5,
color= "grey40",
inherit.aes = FALSE
) +
geom_text(
data = label_df,
aes(
x = as.numeric(factor(trt)) -1.3* dodge_amt,
y = shoot_lab_y,
label = shoot_cld
),
inherit.aes = FALSE,
size = 5
) +
geom_text(
data = label_df,
aes(
x = as.numeric(factor(trt)) + 1.3*dodge_amt,
y = root_lab_y,
label = toupper(root_cld)
),
inherit.aes = FALSE,
size = 3.8,
color= "grey40"
) +
facet_grid(vars(stage), vars(location), scales = "free_y") +
labs(
x = "Fertilizer-N treatment",
y = expression("N derived from fertilizer in crop biomass (kg ha"^{-1}*")"),
fill = "Component"
) +
scale_fill_manual(
values = c(
"Root" = "#D8A24A",
"Shoot" = "#4A9B83"
)
)+
theme_Manjot() +
theme(axis.text.x = element_text(angle = 45, hjust = 1))
p1
ggsave("Ndffbiomass3.png", p1, device = ragg::agg_png, width = 12, height = 8, units = "in", dpi = 300)15N recovery % in all pools (except whole soil) graph
plotdata <- InF %>%
dplyr::select(
trt,
location,
stage,
row,
per_Recov_MAOM_tot,
per_Recov_TDN_tot,
per_Recov_MBN_tot,
perRecov_shoot,
perRecov_root
) %>%
dplyr::group_by(trt, location, stage) %>%
dplyr::summarise(
mean_arootrec = mean(perRecov_root, na.rm = TRUE),
mean_bshootrec = mean(perRecov_shoot, na.rm = TRUE),
mean_maomrec = mean(per_Recov_MAOM_tot, na.rm = TRUE),
mean_mbnrec = mean(
per_Recov_MBN_tot,
na.rm = TRUE
),
mean_tdnrec = mean(per_Recov_TDN_tot, na.rm = TRUE),
.groups = "drop"
) %>%
tidyr::pivot_longer(
cols = c(
mean_arootrec,
mean_bshootrec,
mean_maomrec,
mean_mbnrec,
mean_tdnrec
),
names_to = "pool",
values_to = "amount"
) %>%
dplyr::mutate(
pool = factor(
pool,
levels = c(
"mean_arootrec",
"mean_bshootrec",
"mean_maomrec",
"mean_mbnrec",
"mean_tdnrec"
)
)
)
# Invisible values used to set separate y-axis ranges for V3 and V6
axis_caps <- plotdata %>%
dplyr::distinct(stage, location) %>%
dplyr::mutate(
trt = dplyr::first(plotdata$trt),
amount = dplyr::if_else(stage == "V3", 110, 60)
)
p1 <- ggplot2::ggplot(
plotdata,
ggplot2::aes(
x = factor(trt),
y = amount,
fill = pool
)
) +
ggplot2::geom_col(
position = "stack",
width = 0.60
) +
# Invisible points force V3 to 110 and V6 to 60
ggplot2::geom_blank(
data = axis_caps,
ggplot2::aes(
x = factor(trt),
y = amount
),
inherit.aes = FALSE
) +
ggplot2::facet_grid(
ggplot2::vars(stage),
ggplot2::vars(location),
scales = "free_y"
) +
ggplot2::scale_fill_manual(
values = c(
"mean_arootrec" = "#D8A24A", # Root: muted gold
"mean_bshootrec" = "#4A9B83", # Shoot: muted green
"mean_maomrec" = "#8E5A8A", # MAOM: warm purple
"mean_mbnrec" = "#C48FB3", # MBN: light warm purple
"mean_tdnrec" = "#6678B8" # TDN: cool blue-purple
),
breaks = c(
"mean_arootrec",
"mean_bshootrec",
"mean_maomrec",
"mean_mbnrec",
"mean_tdnrec"
),
labels = c(
"Root",
"Shoot",
"MAOM",
"MBN",
"TDN"
)
) +
ggplot2::labs(
x = "Fertilizer treatment",
y = expression({}^{15}*"N recovery (%)"),
fill = "Pools"
) +
theme_Manjot() +
ggplot2::theme(
axis.text = ggplot2::element_text(size = 15),
axis.text.x = ggplot2::element_text(
angle = 45,
hjust = 1
),
axis.title = ggplot2::element_text(size = 16),
plot.title = ggplot2::element_text(
size = 18,
face = "bold"
),
plot.subtitle = ggplot2::element_text(size = 14),
legend.title = ggplot2::element_text(size = 15),
legend.text = ggplot2::element_text(size = 14)
)
p1
ggsave("15Nrecoveryall.png", p1, device = ragg::agg_png, width = 12, height = 8, units = "in", dpi = 300)Band root biomass graph
model_v6_neb_root <- lm(root_wt ~ trt + row, data = dat_Neb_V6_Band)
model_v6_nh_root <- lm(root_wt ~ trt + row, data = dat_NH_V6_Band)
emm_v6_neb_root <- emmeans(model_v6_neb_root, ~ trt)
emm_v6_nh_root <- emmeans(model_v6_nh_root, ~ trt)
cld_v6_neb_root <- cld(emm_v6_neb_root, Letters = letters, alpha = 0.05, reversed= TRUE)
cld_v6_nh_root <- cld(emm_v6_nh_root, Letters = letters, alpha = 0.05, reversed= TRUE)
root_cld_df <- bind_rows(
as.data.frame(cld_v6_neb_root) %>% mutate(stage = "V6", location = "NE"),
as.data.frame(cld_v6_nh_root) %>% mutate(stage = "V6", location = "NH")
) %>%
dplyr:: rename(root_cld = .group) %>%
dplyr:: mutate(root_cld = gsub(" ", "", root_cld)) %>%
dplyr:: select(trt, location, root_cld)
plot_summary <- Band_V6 %>%
dplyr::select(trt, location, row, root_wt ) %>%
dplyr::group_by(trt, location) %>%
dplyr:: summarise(
mean_root = mean(root_wt, na.rm = TRUE),
se_root = sd(root_wt, na.rm = TRUE) / sqrt(sum(!is.na(root_wt ))),
.groups = "drop"
) %>%
left_join(root_cld_df, by = c("trt", "location"))
plotdata <- plot_summary %>%
dplyr:: select(trt, location, mean_root) %>%
pivot_longer(
cols = c(mean_root),
names_to = "pool",
values_to = "amount"
) %>%
dplyr:: mutate(
pool = factor(pool,
levels = c("mean_root"),
labels = c("Root"))
)
# root error bars: centered on top of shoot segment
root_err <- plot_summary %>%
dplyr::mutate(
ymin = mean_root - se_root,
ymax = mean_root + se_root
)
label_df <- plot_summary %>%
dplyr::group_by(location) %>%
dplyr::mutate(
panel_max = max(mean_root + se_root, na.rm = TRUE),
root_lab_y = mean_root + se_root + 0.06 * panel_max
) %>%
ungroup()
dodge_amt <- 0.09
p1 <- ggplot(plotdata, aes(x = factor(trt), y = amount, fill = pool)) +
geom_col(width = 0.60) +
geom_errorbar(
data = root_err,
aes(
x = as.numeric(factor(trt)),
ymin = ymin,
ymax = ymax
),
width = 0.12,
linewidth = 0.5,
inherit.aes = FALSE
) +
geom_text(
data = label_df,
aes(
x = as.numeric(factor(trt)),
y = root_lab_y,
label = root_cld
),
inherit.aes = FALSE,
size = 5
) +
facet_grid(vars(location), scales = "free_y") +
labs(
x = "Fertilizer-N treatment",
y = expression("Band root biomass (kg ha"^{-1}*")"),
fill = "Component"
) +
#scale_fill_brewer(palette = "Pastel2", labels = c("Root")) +
scale_fill_manual(
values = c(
"Root" = "#D8A24A"
)
)+
theme_Manjot() +
theme(axis.text.x = element_text(angle = 45, hjust = 1))
p1
Ndff in TDN graph
#TDN band models
model_v3_neb_shoot <- lm(kgN_dff_TDN_ha_band ~ trt + row, data = dat_Neb_V3)
model_v6_neb_shoot <- lm(kgN_dff_TDN_ha_band~ trt + row, data = dat_Neb_V6)
model_v3_nh_shoot <- lm(kgN_dff_TDN_ha_band~ trt + row, data = dat_NH_V3)
model_v6_nh_shoot <- lm(kgN_dff_TDN_ha_band~ trt + row, data = dat_NH_V6)
emm_v3_neb_shoot <- emmeans(model_v3_neb_shoot, ~ trt)
emm_v6_neb_shoot <- emmeans(model_v6_neb_shoot, ~ trt)
emm_v3_nh_shoot <- emmeans(model_v3_nh_shoot, ~ trt)
emm_v6_nh_shoot <- emmeans(model_v6_nh_shoot, ~ trt)
cld_v3_neb_shoot <- cld(emm_v3_neb_shoot, Letters = letters, alpha = 0.05, sort= TRUE, reversed= TRUE)
cld_v6_neb_shoot <- cld(emm_v6_neb_shoot, Letters = letters, alpha = 0.05, sort= TRUE, reversed= TRUE)
cld_v3_nh_shoot <- cld(emm_v3_nh_shoot, Letters = letters, alpha = 0.05, sort= TRUE, reversed= TRUE)
cld_v6_nh_shoot <- cld(emm_v6_nh_shoot, Letters = letters, alpha = 0.05, sort= TRUE, reversed= TRUE)
shoot_cld_df <- bind_rows(
as.data.frame(cld_v3_neb_shoot) %>% mutate(stage = "V3", location = "NE"),
as.data.frame(cld_v6_neb_shoot) %>% mutate(stage = "V6", location = "NE"),
as.data.frame(cld_v3_nh_shoot) %>% mutate(stage = "V3", location = "NH"),
as.data.frame(cld_v6_nh_shoot) %>% mutate(stage = "V6", location = "NH")
) %>%
dplyr::rename(shoot_cld = .group) %>%
dplyr::mutate(shoot_cld = gsub(" ", "", shoot_cld)) %>%
dplyr::select(trt, location, stage, shoot_cld)
#-----------------------------
# TDN furrow models
#-----------------------------
model_v3_neb_root <- lm(kgN_dff_TDN_ha ~ trt + row, data = dat_Neb_V3_Inf)
model_v6_neb_root <- lm(kgN_dff_TDN_ha~ trt + row, data = dat_Neb_V6_Inf)
model_v3_nh_root <- lm(kgN_dff_TDN_ha ~ trt + row, data = dat_NH_V3_Inf)
model_v6_nh_root <- lm(kgN_dff_TDN_ha ~ trt + row, data = dat_NH_V6_Inf)
emm_v3_neb_root <- emmeans(model_v3_neb_root, ~ trt)
emm_v6_neb_root <- emmeans(model_v6_neb_root, ~ trt)
emm_v3_nh_root <- emmeans(model_v3_nh_root, ~ trt)
emm_v6_nh_root <- emmeans(model_v6_nh_root, ~ trt)
cld_v3_neb_root <- cld(emm_v3_neb_root, Letters = letters, alpha = 0.05, sort= TRUE, reversed= TRUE)
cld_v6_neb_root <- cld(emm_v6_neb_root, Letters = letters, alpha = 0.05, sort= TRUE, reversed= TRUE)
cld_v3_nh_root <- cld(emm_v3_nh_root, Letters = letters, alpha = 0.05, sort= TRUE, reversed= TRUE)
cld_v6_nh_root <- cld(emm_v6_nh_root, Letters = letters, alpha = 0.05, sort= TRUE, reversed= TRUE)
root_cld_df <- bind_rows(
as.data.frame(cld_v3_neb_root) %>% mutate(stage = "V3", location = "NE"),
as.data.frame(cld_v6_neb_root) %>% mutate(stage = "V6", location = "NE"),
as.data.frame(cld_v3_nh_root) %>% mutate(stage = "V3", location = "NH"),
as.data.frame(cld_v6_nh_root) %>% mutate(stage = "V6", location = "NH")
) %>%
dplyr:: rename(root_cld = .group) %>%
dplyr:: mutate(root_cld = gsub(" ", "", root_cld)) %>%
dplyr:: select(trt, location, stage, root_cld)
plot_summary <- InF %>%
dplyr::select(trt, location, stage, row, kgN_dff_TDN_ha_band, kgN_dff_TDN_ha) %>%
dplyr::group_by(trt, location, stage) %>%
dplyr:: summarise(
mean_root = mean(kgN_dff_TDN_ha, na.rm = TRUE),
se_root = sd(kgN_dff_TDN_ha, na.rm = TRUE) / sqrt(sum(!is.na(kgN_dff_TDN_ha))),
mean_shoot = mean(kgN_dff_TDN_ha_band, na.rm = TRUE),
se_shoot = sd(kgN_dff_TDN_ha_band, na.rm = TRUE) / sqrt(sum(!is.na(kgN_dff_TDN_ha_band))),
.groups = "drop"
) %>%
left_join(shoot_cld_df, by = c("trt", "location", "stage")) %>%
left_join(root_cld_df, by = c("trt", "location", "stage"))
plotdata <- plot_summary %>%
dplyr:: select(trt, location, stage, mean_root, mean_shoot) %>%
pivot_longer(
cols = c(mean_root, mean_shoot),
names_to = "pool",
values_to = "amount"
) %>%
dplyr:: mutate(
pool = factor(pool,
levels = c("mean_root", "mean_shoot"),
labels = c("Furrow", "Band"))
)
#=========================================================
# Prepare raw TDN data in long format
#=========================================================
trt_levels <- c(
"InF-5",
"InF-10",
"UAN-74",
"Both-79",
"UAN-148",
"Both-158"
)
tdn_plot_data <- InF %>%
dplyr::select(
trt,
location,
stage,
row,
kgN_dff_TDN_ha_band,
kgN_dff_TDN_ha
) %>%
tidyr::pivot_longer(
cols = c(
kgN_dff_TDN_ha_band,
kgN_dff_TDN_ha
),
names_to = "position_plot",
values_to = "amount"
) %>%
dplyr::mutate(
position_plot = dplyr::recode(
position_plot,
"kgN_dff_TDN_ha_band" = "Band",
"kgN_dff_TDN_ha" = "Furrow"
),
trt = factor(
as.character(trt),
levels = trt_levels
),
stage = factor(
as.character(stage),
levels = c("V3", "V6")
),
location = factor(
as.character(location),
levels = c("NE", "NH")
),
position_plot = factor(
position_plot,
levels = c("Band", "Furrow")
)
) %>%
dplyr::filter(
stage %in% c("V3", "V6"),
location %in% c("NE", "NH"),
!is.na(amount)
)
#=========================================================
# Combine CLD letters for Band and Furrow
#=========================================================
position_cld_df <- dplyr::bind_rows(
shoot_cld_df %>%
dplyr::transmute(
trt = as.character(trt),
location = as.character(location),
stage = as.character(stage),
position_plot = "Band",
cld = shoot_cld
),
root_cld_df %>%
dplyr::transmute(
trt = as.character(trt),
location = as.character(location),
stage = as.character(stage),
position_plot = "Furrow",
cld = root_cld
)
) %>%
dplyr::mutate(
trt = factor(trt, levels = trt_levels),
location = factor(
location,
levels = c("NE", "NH")
),
stage = factor(
stage,
levels = c("V3", "V6")
),
position_plot = factor(
position_plot,
levels = c("Band", "Furrow")
)
)
#=========================================================
# Calculate boxplot upper whiskers for CLD positions
#=========================================================
label_df <- tdn_plot_data %>%
dplyr::group_by(
location,
position_plot,
stage,
trt
) %>%
dplyr::summarise(
y_pos = boxplot.stats(amount)$stats[5],
.groups = "drop"
) %>%
dplyr::left_join(
position_cld_df,
by = c(
"trt",
"location",
"stage",
"position_plot"
)
) %>%
dplyr::group_by(
location,
position_plot
) %>%
dplyr::mutate(
panel_min = min(y_pos, na.rm = TRUE),
panel_max = max(y_pos, na.rm = TRUE),
panel_range = panel_max - panel_min,
label_offset = dplyr::if_else(
panel_range > 0,
0.04 * panel_range,
0.04 * abs(panel_max)
),
y_pos = y_pos + label_offset
) %>%
dplyr::ungroup()
#=========================================================
# Create dodged boxplot
#=========================================================
dodge_width <- 0.75
tdn_plot_data_plot <- tdn_plot_data %>%
dplyr::group_by(
location,
position_plot,
stage,
trt
) %>%
dplyr::ungroup()
label_df <- tdn_plot_data_plot %>%
dplyr::group_by(
location,
position_plot,
stage,
trt
) %>%
dplyr::summarise(
y_pos = boxplot.stats(amount)$stats[5],
.groups = "drop"
) %>%
dplyr::left_join(
position_cld_df,
by = c(
"trt",
"location",
"stage",
"position_plot"
)
) %>%
dplyr::group_by(location, position_plot) %>%
dplyr::mutate(
panel_max = max(y_pos, na.rm = TRUE),
label_offset = 0.05 * pmax(panel_max, 0.1),
y_pos = y_pos + label_offset
) %>%
dplyr::ungroup()
p1 <- ggplot2::ggplot(
tdn_plot_data_plot,
ggplot2::aes(
x = trt,
y = amount,
fill = stage,
group = interaction(trt, stage)
)
) +
ggplot2::geom_boxplot(
position = ggplot2::position_dodge(
width = dodge_width
),
width = 0.60,
linewidth = 0.5
) +
ggplot2::geom_text(
data = label_df,
ggplot2::aes(
x = trt,
y = y_pos,
label = cld,
group = stage
),
position = ggplot2::position_dodge(
width = dodge_width
),
inherit.aes = FALSE,
size = 5
) +
# Band and Furrow receive separate y-axis scales
# NE and NH remain as columns
ggplot2::facet_grid(
rows = ggplot2::vars(position_plot),
cols = ggplot2::vars(location),
scales = "free_y"
) +
# Blue palette consistent with TDN
ggplot2::scale_fill_manual(
values = c(
"V3" = "#4F66A8",
"V6" = "#8DBCE6"
),
breaks = c("V3", "V6")
) +
ggplot2::scale_y_continuous(
expand = ggplot2::expansion(
mult = c(0.03, 0.15)
)
) +
ggplot2::labs(
x = "Fertilizer-N treatment",
y = expression(
"Amount N dff in TDN (kg N ha"^{-1} * ")"
),
fill = "Growth stage"
) +
theme_Manjot() +
ggplot2::theme(
axis.text = ggplot2::element_text(size = 13),
axis.text.x = ggplot2::element_text(
angle = 45,
hjust = 1
),
axis.title = ggplot2::element_text(size = 16),
plot.title = ggplot2::element_text(
size = 18,
face = "bold"
),
legend.title = ggplot2::element_text(size = 14),
legend.text = ggplot2::element_text(size = 13),
strip.text = ggplot2::element_text(
size = 14,
face = "bold"
)
)
p1
ggsave("TDNexp2dffbox.jpg", p1, device = ragg::agg_png, width = 12, height = 8, units = "in", dpi = 300)DIN graph
#-------------
# DIN band models
#--------------
model_v3_neb_shoot <- lm(ug_DIN_g~ trt + row, data = dat_Neb_V3_Band)
model_v6_neb_shoot <- lm(ug_DIN_g~ trt + row, data = dat_Neb_V6_Band)
model_v3_nh_shoot <- lm(ug_DIN_g~ trt + row, data = dat_NH_V3_Band)
model_v6_nh_shoot <- lm(ug_DIN_g~ trt + row, data = dat_NH_V6_Band)
emm_v3_neb_shoot <- emmeans(model_v3_neb_shoot, ~ trt)
emm_v6_neb_shoot <- emmeans(model_v6_neb_shoot, ~ trt)
emm_v3_nh_shoot <- emmeans(model_v3_nh_shoot, ~ trt)
emm_v6_nh_shoot <- emmeans(model_v6_nh_shoot, ~ trt)
cld_v3_neb_shoot <- cld(emm_v3_neb_shoot, Letters = letters, alpha = 0.05, sort= TRUE, reversed= TRUE)
cld_v6_neb_shoot <- cld(emm_v6_neb_shoot, Letters = letters, alpha = 0.05, sort= TRUE, reversed= TRUE)
cld_v3_nh_shoot <- cld(emm_v3_nh_shoot, Letters = letters, alpha = 0.05, sort= TRUE, reversed= TRUE)
cld_v6_nh_shoot <- cld(emm_v6_nh_shoot, Letters = letters, alpha = 0.05, sort= TRUE, reversed= TRUE)
shoot_cld_df <- bind_rows(
as.data.frame(cld_v3_neb_shoot) %>% mutate(stage = "V3", location = "NE"),
as.data.frame(cld_v6_neb_shoot) %>% mutate(stage = "V6", location = "NE"),
as.data.frame(cld_v3_nh_shoot) %>% mutate(stage = "V3", location = "NH"),
as.data.frame(cld_v6_nh_shoot) %>% mutate(stage = "V6", location = "NH")
) %>%
dplyr::rename(shoot_cld = .group) %>%
dplyr::mutate(shoot_cld = gsub(" ", "", shoot_cld)) %>%
dplyr::select(trt, location, stage, shoot_cld)
#-----------------------------
# DIN Furrow models
#-----------------------------
model_v3_neb_root <- lm(ug_DIN_g ~ trt + row, data = dat_Neb_V3_Inf)
model_v6_neb_root <- lm(ug_DIN_g~ trt + row, data = dat_Neb_V6_Inf)
model_v3_nh_root <- lm(ug_DIN_g ~ trt + row, data = dat_NH_V3_Inf)
model_v6_nh_root <- lm(ug_DIN_g ~ trt + row, data = dat_NH_V6_Inf)
emm_v3_neb_root <- emmeans(model_v3_neb_root, ~ trt)
emm_v6_neb_root <- emmeans(model_v6_neb_root, ~ trt)
emm_v3_nh_root <- emmeans(model_v3_nh_root, ~ trt)
emm_v6_nh_root <- emmeans(model_v6_nh_root, ~ trt)
cld_v3_neb_root <- cld(emm_v3_neb_root, Letters = letters, alpha = 0.05, sort= TRUE, reversed= TRUE)
cld_v6_neb_root <- cld(emm_v6_neb_root, Letters = letters, alpha = 0.05, sort= TRUE, reversed= TRUE)
cld_v3_nh_root <- cld(emm_v3_nh_root, Letters = letters, alpha = 0.05, sort= TRUE, reversed= TRUE)
cld_v6_nh_root <- cld(emm_v6_nh_root, Letters = letters, alpha = 0.05, sort= TRUE, reversed= TRUE)
root_cld_df <- bind_rows(
as.data.frame(cld_v3_neb_root) %>% mutate(stage = "V3", location = "NE"),
as.data.frame(cld_v6_neb_root) %>% mutate(stage = "V6", location = "NE"),
as.data.frame(cld_v3_nh_root) %>% mutate(stage = "V3", location = "NH"),
as.data.frame(cld_v6_nh_root) %>% mutate(stage = "V6", location = "NH")
) %>%
dplyr:: rename(root_cld = .group) %>%
dplyr:: mutate(root_cld = gsub(" ", "", root_cld)) %>%
dplyr:: select(trt, location, stage, root_cld)
position_cld_df <- dplyr::bind_rows(
shoot_cld_df %>%
dplyr::transmute(
trt,
location,
stage,
position = "Band",
cld = shoot_cld
),
root_cld_df %>%
dplyr::transmute(
trt,
location,
stage,
position = "In-row",
cld = root_cld
)
)
din_plot_data <- data_new %>%
dplyr::filter(stage %in% c("V3", "V6")) %>%
dplyr::mutate(
position_plot = dplyr::case_when(
position == "Band" ~ "Band",
position %in% c("Inf", "InF", "Furrow", "In-row") ~ "Furrow",
TRUE ~ NA_character_
)
) %>%
dplyr::filter(
!is.na(position_plot),
!is.na(ug_DIN_g)
) %>%
dplyr::mutate(
position_plot = factor(
position_plot,
levels = c("Band", "Furrow")
)
)
#---------------------------------------------------------
# Combine CLD letters for the two positions
#---------------------------------------------------------
position_cld_df <- dplyr::bind_rows(
shoot_cld_df %>%
dplyr::transmute(
trt,
location,
stage,
position_plot = "Band",
cld = shoot_cld
),
root_cld_df %>%
dplyr::transmute(
trt,
location,
stage,
position_plot = "Furrow",
cld = root_cld
)
) %>%
dplyr::mutate(
position_plot = factor(
position_plot,
levels = c("Band", "Furrow")
)
)
#---------------------------------------------------------
# Calculate upper-whisker positions for CLD letters
#---------------------------------------------------------
label_df <- din_plot_data %>%
dplyr::group_by(
location,
position_plot,
stage,
trt
) %>%
dplyr::summarise(
y_pos = boxplot.stats(ug_DIN_g)$stats[5],
.groups = "drop"
) %>%
dplyr::left_join(
position_cld_df,
by = c(
"trt",
"location",
"stage",
"position_plot"
)
) %>%
dplyr::group_by(location, position_plot, stage) %>%
dplyr::mutate(
stage_min = min(y_pos, na.rm = TRUE),
stage_max = max(y_pos, na.rm = TRUE),
stage_range = stage_max - stage_min,
label_offset = 0.08 * pmax(stage_range, stage_max, 1),
y_pos = y_pos + label_offset
) %>%
dplyr::ungroup()
#---------------------------------------------------------
# Dodged DIN boxplot
#---------------------------------------------------------
dodge_width <- 0.75
p1 <- ggplot2::ggplot(
din_plot_data,
ggplot2::aes(
x = factor(trt),
y = ug_DIN_g,
fill = stage,
group = interaction(trt, stage)
)
) +
ggplot2::geom_boxplot(
position = ggplot2::position_dodge(width = dodge_width),
width = 0.60,
linewidth = 0.5
) +
ggplot2::geom_text(
data = label_df,
ggplot2::aes(
x = factor(trt),
y = y_pos,
label = cld,
group = stage
),
position = ggplot2::position_dodge(width = dodge_width),
inherit.aes = FALSE,
size = 5
) +
# Band and Furrow are rows; NE and NH are columns
ggplot2::facet_grid(
rows = ggplot2::vars(position_plot),
cols = ggplot2::vars(location),
scales = "free_y",
labeller = ggplot2::labeller(
position_plot = c(
Band = "Band",
Furrow = "Furrow"
)
)
) +
ggplot2::scale_fill_manual(
values = c(
"V3" = "#4F66A8",
"V6" = "#8DBCE6"
)
) +
ggplot2::scale_y_continuous(
expand = ggplot2::expansion(mult = c(0.03, 0.15))
) +
ggplot2::labs(
x = "Fertilizer-N treatment",
y = expression(
"DIN (" * mu * "g N g"^{-1} * " soil)"
),
fill = "Growth stage",
) +
theme_Manjot() +
ggplot2::theme(
axis.text = ggplot2::element_text(size = 13),
axis.text.x = ggplot2::element_text(
angle = 45,
hjust = 1
),
axis.title = ggplot2::element_text(size = 16),
legend.title = ggplot2::element_text(size = 14),
legend.text = ggplot2::element_text(size = 13),
strip.text = ggplot2::element_text(
size = 14,
face = "bold"
)
)
p1