Prepare time-series data
# Create long-format data with all time points
dw_raw <- pheno %>%
mutate(STW121 = STWHV) %>%
select(
plotId,
Genotype,
Genotype_label,
Treatment,
starts_with("STW")
) %>%
pivot_longer(
cols = c("STW40", "STW50", "STW60", "STW121"),
names_to = "DAP_string",
values_to = "STW"
) %>%
mutate(
DAP = as.integer(gsub("STW", "", DAP_string)),
# Ensure factor levels are maintained
Genotype = factor(Genotype, levels = c("CTRL", "Inv4m")),
Genotype_label = factor(Genotype_label, levels = c("CTRL", "<i>Inv4m</i>"))
) %>%
filter(!is.na(STW)) %>%
droplevels()
cat("Time-series data dimensions:", dim(dw_raw), "\n")
## Time-series data dimensions: 236 8
cat("Time points:", sort(unique(dw_raw$DAP)), "\n")
## Time points: 40 50 60 121
cat("Genotype_label levels:", levels(dw_raw$Genotype_label), "\n")
## Genotype_label levels: CTRL <i>Inv4m</i>
Time-point specific t-tests
# Perform t-tests at each time point within each genotype
treatment_timeseries_test <- dw_raw %>%
ungroup() %>%
filter(DAP > 1) %>%
mutate(DAP = as.factor(DAP)) %>%
group_by(Genotype, DAP) %>%
t_test(STW ~ Treatment, ref.group = "+P") %>%
adjust_pvalue(method = "fdr") %>%
add_significance() %>%
add_y_position(scales = "free_y") %>%
add_x_position(x = "DAP", group = NULL, dodge = 1) %>%
mutate(
Genotype_label = factor(
case_when(
Genotype == "CTRL" ~ "CTRL",
Genotype == "Inv4m" ~ "<i>Inv4m</i>",
TRUE ~ as.character(Genotype)
),
levels = c("CTRL", "<i>Inv4m</i>")
)
)
# Adjust bracket positions relative to data range (not absolute units)
data_range <- range(dw_raw$STW[dw_raw$DAP > 1], na.rm = TRUE)
y_range <- diff(data_range)
treatment_timeseries_test <- treatment_timeseries_test %>%
mutate(y.position = y.position + 0.02 * y_range) # 8% of data range
cat("\n=== Time-series t-test summary ===\n")
##
## === Time-series t-test summary ===
print(treatment_timeseries_test %>%
select(Genotype, DAP, statistic, p, p.adj, p.adj.signif))
## # A tibble: 8 × 6
## Genotype DAP statistic p p.adj p.adj.signif
## <fct> <fct> <dbl> <dbl> <dbl> <chr>
## 1 CTRL 40 5.93 0.0000039 0.0000078 ****
## 2 CTRL 50 10.2 0.0000000324 0.000000130 ****
## 3 CTRL 60 7.18 0.000000961 0.00000256 ****
## 4 CTRL 121 4.73 0.0000995 0.000159 ***
## 5 Inv4m 40 4.01 0.000445 0.000593 ***
## 6 Inv4m 50 9.49 0.0000000271 0.000000130 ****
## 7 Inv4m 60 3.97 0.00141 0.00161 **
## 8 Inv4m 121 3.25 0.00369 0.00369 **
Visualize stover weight trajectory
pd <- position_dodge(1)
p_dw <- dw_raw %>%
ungroup() %>%
filter(DAP > 1) %>%
ggplot(aes(x = as.factor(DAP), y = STW, col = Treatment)) +
ggtitle("Stover Dry Weight") +
ylab("g") +
xlab("Days After Planting") +
geom_boxplot(
width = 0.25,
linewidth = 1,
alpha = 0,
position = pd
) %>%
with_shadow(
colour = "black",
x_offset = 0,
y_offset = 0,
sigma = 1
) +
geom_quasirandom(
size = 2,
width = 0.25,
dodge.width = 1
) %>%
with_shadow(
colour = "black",
x_offset = 0,
y_offset = 0,
sigma = 1
) +
facet_wrap(. ~ Genotype_label) +
stat_pvalue_manual(
treatment_timeseries_test,
size = 10,
bracket.size = 0.8,
tip.length = 0.01,
hide.ns = TRUE
) +
scale_color_manual(values = pal) +
scale_y_continuous(breaks = seq(0,150,25),
expand = expansion(mult = c(0.05, 0.1))) +
pheno_theme3 +
theme(
legend.position = c(0.12, 0.7),
strip.text = element_markdown()
)
print(p_dw)
# Extract Inv4m genotype effects for all stover weight traits with CIs
# Extract both Inv4m genotype effects and -P treatment effects for stover weights
stover_effects <- effects %>%
filter(
predictor %in% c("GenotypeInv4m", "Treatment-P"),
grepl("^STW", response)
) %>%
# Calculate reference means
left_join(
effects %>%
filter(grepl("Intercept|Treatment-P|GenotypeInv4m", predictor)) %>%
select(response, predictor, Estimate) %>%
pivot_wider(names_from = predictor, values_from = Estimate) %>%
rename(
intercept = `(Intercept)`,
treatment_effect = `Treatment-P`,
genotype_effect = GenotypeInv4m
) %>%
mutate(
# +P marginal mean for treatment effects
plus_p_marginal_mean = intercept + (genotype_effect / 2),
# CTRL marginal mean for genotype effects
ctrl_marginal_mean = intercept + (treatment_effect / 2)
) %>%
select(response, plus_p_marginal_mean, ctrl_marginal_mean),
by = "response"
) %>%
mutate(
# Use appropriate denominator based on effect type
reference_mean = if_else(
predictor == "Treatment-P",
plus_p_marginal_mean,
ctrl_marginal_mean
),
# Calculate 95% CI
ci_lower = Estimate - 1.96 * Std.Error,
ci_upper = Estimate + 1.96 * Std.Error,
# Calculate percentage and its CI
estimate_pct = (Estimate / reference_mean) * 100,
pct_ci_lower = (ci_lower / reference_mean) * 100,
pct_ci_upper = (ci_upper / reference_mean) * 100,
# Extract time point
DAP = case_when(
response == "STW40" ~ "40",
response == "STW50" ~ "50",
response == "STW60" ~ "60",
response == "STWHV" ~ "121",
TRUE ~ NA_character_
),
DAP = factor(DAP, levels = c("40", "50", "60", "121")),
# Clean predictor names
effect_type = case_when(
predictor == "GenotypeInv4m" ~ "Inv4m",
predictor == "Treatment-P" ~ "-P",
TRUE ~ predictor
),
effect_type = factor(effect_type, levels = c("Inv4m", "-P")),
significant = p.adjust < 0.05
) %>%
arrange(DAP, effect_type)
# Fit regression for -P data
minus_p_data <- stover_effects %>%
filter(effect_type == "-P") %>%
mutate(DAP_numeric = as.numeric(as.character(DAP)))
lm_minus_p <- lm(estimate_pct ~ DAP_numeric, data = minus_p_data)
r_squared_minus_p <- summary(lm_minus_p)$r.squared
p_value_minus_p <- summary(lm_minus_p)$coefficients[2, 4]
# Fit regression for Inv4m data
inv4m_data <- stover_effects %>%
filter(effect_type == "Inv4m") %>%
mutate(DAP_numeric = as.numeric(as.character(DAP)))
lm_inv4m <- lm(estimate_pct ~ DAP_numeric, data = inv4m_data)
r_squared_inv4m <- summary(lm_inv4m)$r.squared
p_value_inv4m <- summary(lm_inv4m)$coefficients[2, 4]
# Create annotation text for both
regression_label_minus_p <- sprintf(
"italic(R)^2 == %.3f~~italic(p) == %.3f",
r_squared_minus_p,
p_value_minus_p
)
regression_label_inv4m <- sprintf(
"italic(R)^2 == %.3f~~italic(p) == %.3f",
r_squared_inv4m,
p_value_inv4m
)
# Create plot comparing both effects
p_biomass_comparison <- stover_effects %>%
mutate(DAP_numeric = as.numeric(as.character(DAP))) %>%
ggplot(aes(x = DAP_numeric, y = -estimate_pct, group = effect_type)) +
# Add regression line for -P
geom_smooth(
data = minus_p_data,
method = "lm",
se = FALSE,
linewidth = 1,
linetype = "dashed",
color = "black"
) +
# Add regression line for Inv4m
geom_smooth(
data = inv4m_data,
method = "lm",
se = FALSE,
linewidth = 1,
linetype = "dotted",
color = "black"
) +
geom_pointrange(
aes(ymin = -pct_ci_lower, ymax = -pct_ci_upper,
shape = effect_type),
color = "black",
size = 1,
linewidth = 1,
fatten = 5,
position = position_dodge(width = 3)
) +
# Add R² and p-value annotation for -P
annotate(
"text",
x = 80,
y = 45,
label = regression_label_minus_p,
parse = TRUE,
size = 6,
hjust = 0
) +
# Add R² and p-value annotation for Inv4m
annotate(
"text",
x = 80,
y = 5,
label = regression_label_inv4m,
parse = TRUE,
size = 6,
hjust = 0
) +
scale_x_continuous(
breaks = c(40, 50, 60, 121),
labels = c("40", "50", "60", "121")
) +
scale_shape_manual(
values = c(18, 16),
labels = c("<i>Inv4m</i>", "-P")
) +
guides(shape = guide_legend(reverse = TRUE)) +
labs(
title = "Stover Biomass Reduction",
x = "Days After Planting",
y = "% reference mean",
shape = "Predictor"
) +
theme_classic2(base_size = 20) +
theme(
plot.title = element_markdown(hjust = 0.5, face = "bold"),
axis.text.x = element_text(color = "black", size = 18),
legend.position = c(0.8, 0.8),
legend.text = element_markdown(size = 20)
)
print(p_biomass_comparison)
p_stw <- ggarrange(
p_dw,
p_biomass_comparison,
nrow = 1,
ncol = 2,
labels = c("E", "F"),
font.label = list(size = 30)
)
p_stw
