##Import data and clean
library(dplyr)
##
## Attaching package: 'dplyr'
## The following objects are masked from 'package:stats':
##
## filter, lag
## The following objects are masked from 'package:base':
##
## intersect, setdiff, setequal, union
library(tidyr)
library(ggplot2)
library(scales)
data <- read.csv(
"exp04_cut_soybeans_sc1_20260729.csv",
check.names = FALSE
)
# Inspect the original dataset
names(data)
## [1] "Date" "Plant_ID"
## [3] "Cut_Status" "Water_Status"
## [5] "Lowest_Mature_Meristem_SC" "Highest_Mature_Meristem_SC"
## [7] "Notes"
head(data)
## Date Plant_ID Cut_Status Water_Status Lowest_Mature_Meristem_SC
## 1 7/29/2026 1a Cut Watered 198.1
## 2 7/29/2026 1b Intact Watered 64.9
## 3 7/29/2026 2a Cut Drought 168.5
## 4 7/29/2026 2b Intact Drought 49.2
## 5 7/29/2026 3a Cut Watered 273.4
## 6 7/29/2026 3b Intact Watered 188.7
## Highest_Mature_Meristem_SC Notes
## 1 220.9
## 2 215.4
## 3 130.8
## 4 110.0
## 5 305.8
## 6 286.4
str(data)
## 'data.frame': 42 obs. of 7 variables:
## $ Date : chr "7/29/2026" "7/29/2026" "7/29/2026" "7/29/2026" ...
## $ Plant_ID : chr "1a" "1b" "2a" "2b" ...
## $ Cut_Status : chr "Cut" "Intact" "Cut" "Intact" ...
## $ Water_Status : chr "Watered" "Watered" "Drought" "Drought" ...
## $ Lowest_Mature_Meristem_SC : num 198.1 64.9 168.5 49.2 273.4 ...
## $ Highest_Mature_Meristem_SC: num 221 215 131 110 306 ...
## $ Notes : chr "" "" "" "" ...
data_clean <- data %>%
mutate(
Date = as.Date(
Date,
format = "%m/%d/%Y"
),
Architecture = case_when(
Cut_Status == "Cut" ~ "Compact",
Cut_Status == "Intact" ~ "Erect",
TRUE ~ NA_character_
),
Architecture = factor(
Architecture,
levels = c("Erect", "Compact")
),
Water_Status = factor(
Water_Status,
levels = c("Watered", "Drought")
),
Plant_ID = factor(Plant_ID)
)
# Convert the highest- and lowest-position columns into long format
analysis_data <- data_clean %>%
pivot_longer(
cols = c(
Lowest_Mature_Meristem_SC,
Highest_Mature_Meristem_SC
),
names_to = "Mature_Position",
values_to = "SC"
) %>%
mutate(
Mature_Position = recode(
Mature_Position,
"Lowest_Mature_Meristem_SC" = "Lowest Mature",
"Highest_Mature_Meristem_SC" = "Highest Mature"
),
Mature_Position = factor(
Mature_Position,
levels = c(
"Lowest Mature",
"Highest Mature"
)
),
SC = as.numeric(
gsub("[^0-9.-]", "", SC)
)
)
# Confirm that 4a and 4b remain in the analysis
analysis_data %>%
filter(Plant_ID %in% c("4a", "4b")) %>%
arrange(Plant_ID, Date, Mature_Position)
## # A tibble: 12 × 8
## Date Plant_ID Cut_Status Water_Status Notes Architecture
## <date> <fct> <chr> <fct> <chr> <fct>
## 1 2026-07-29 4a Cut Drought "Mangled leaf noted… Compact
## 2 2026-07-29 4a Cut Drought "Mangled leaf noted… Compact
## 3 2026-07-30 4a Cut Drought "" Compact
## 4 2026-07-30 4a Cut Drought "" Compact
## 5 2026-07-31 4a Cut Drought "" Compact
## 6 2026-07-31 4a Cut Drought "" Compact
## 7 2026-07-29 4b Intact Watered "" Erect
## 8 2026-07-29 4b Intact Watered "" Erect
## 9 2026-07-30 4b Intact Watered "" Erect
## 10 2026-07-30 4b Intact Watered "" Erect
## 11 2026-07-31 4b Intact Watered "" Erect
## 12 2026-07-31 4b Intact Watered "" Erect
## # ℹ 2 more variables: Mature_Position <fct>, SC <dbl>
# Check sample sizes
analysis_data %>%
count(
Date,
Architecture,
Water_Status,
Mature_Position
)
## # A tibble: 24 × 5
## Date Architecture Water_Status Mature_Position n
## <date> <fct> <fct> <fct> <int>
## 1 2026-07-29 Erect Watered Lowest Mature 4
## 2 2026-07-29 Erect Watered Highest Mature 4
## 3 2026-07-29 Erect Drought Lowest Mature 3
## 4 2026-07-29 Erect Drought Highest Mature 3
## 5 2026-07-29 Compact Watered Lowest Mature 3
## 6 2026-07-29 Compact Watered Highest Mature 3
## 7 2026-07-29 Compact Drought Lowest Mature 4
## 8 2026-07-29 Compact Drought Highest Mature 4
## 9 2026-07-30 Erect Watered Lowest Mature 4
## 10 2026-07-30 Erect Watered Highest Mature 4
## # ℹ 14 more rows
# Check for missing values
analysis_data %>%
filter(is.na(SC))
## # A tibble: 0 × 8
## # ℹ 8 variables: Date <date>, Plant_ID <fct>, Cut_Status <chr>,
## # Water_Status <fct>, Notes <chr>, Architecture <fct>, Mature_Position <fct>,
## # SC <dbl>
##Trajectory Analysis
trajectory_summary <- analysis_data %>%
group_by(
Date,
Architecture,
Water_Status,
Mature_Position
) %>%
summarise(
n = sum(!is.na(SC)),
Mean_SC = mean(
SC,
na.rm = TRUE
),
SD_SC = sd(
SC,
na.rm = TRUE
),
SE_SC = SD_SC / sqrt(n),
.groups = "drop"
)
trajectory_summary
## # A tibble: 24 × 8
## Date Architecture Water_Status Mature_Position n Mean_SC SD_SC
## <date> <fct> <fct> <fct> <int> <dbl> <dbl>
## 1 2026-07-29 Erect Watered Lowest Mature 4 145. 97.0
## 2 2026-07-29 Erect Watered Highest Mature 4 255. 99.6
## 3 2026-07-29 Erect Drought Lowest Mature 3 115. 56.9
## 4 2026-07-29 Erect Drought Highest Mature 3 147. 32.6
## 5 2026-07-29 Compact Watered Lowest Mature 3 200. 72.6
## 6 2026-07-29 Compact Watered Highest Mature 3 206. 108.
## 7 2026-07-29 Compact Drought Lowest Mature 4 130. 35.9
## 8 2026-07-29 Compact Drought Highest Mature 4 178. 54.7
## 9 2026-07-30 Erect Watered Lowest Mature 4 80.8 15.4
## 10 2026-07-30 Erect Watered Highest Mature 4 275. 76.1
## # ℹ 14 more rows
## # ℹ 1 more variable: SE_SC <dbl>
trajectory_plot <- ggplot(
trajectory_summary,
aes(
x = Date,
y = Mean_SC,
color = Architecture,
group = Architecture
)
) +
geom_line(
linewidth = 1.1
) +
geom_point(
size = 3
) +
geom_errorbar(
aes(
ymin = Mean_SC - SE_SC,
ymax = Mean_SC + SE_SC
),
width = 0.08,
linewidth = 0.7
) +
facet_grid(
Water_Status ~ Mature_Position
) +
scale_x_date(
date_labels = "%m/%d",
date_breaks = "1 day"
) +
labs(
title = "Stomatal Conductance Trajectories of Compact and Erect Soybeans",
subtitle = "Mean stomatal conductance ± SE across three measurement dates",
x = "Date",
y = expression(
"Stomatal Conductance (mmol " *
m^-2 * " " * s^-1 * ")"
),
color = "Architecture"
) +
theme_classic(
base_size = 12
) +
theme(
plot.title = element_text(
size = 15,
face = "bold",
hjust = 0.5
),
plot.subtitle = element_text(
size = 10,
hjust = 0.5
),
strip.text = element_text(
size = 11,
face = "bold"
),
axis.title = element_text(
size = 11,
face = "bold"
),
axis.text = element_text(
size = 10
),
legend.position = "bottom"
)
trajectory_plot
###Interpretation
This graph shows whether compact and erect plants maintained stomatal conductance differently over time. It also distinguishes watered and droughted plants and the highest and lowest mature positions.
##Drought Penalty
penalty_summary <- analysis_data %>%
group_by(
Date,
Architecture,
Mature_Position,
Water_Status
) %>%
summarise(
n = sum(!is.na(SC)),
Mean_SC = mean(
SC,
na.rm = TRUE
),
SD_SC = sd(
SC,
na.rm = TRUE
),
SE_SC = SD_SC / sqrt(n),
.groups = "drop"
)
drought_penalty_data <- penalty_summary %>%
select(
Date,
Architecture,
Mature_Position,
Water_Status,
Mean_SC,
SE_SC
) %>%
pivot_wider(
names_from = Water_Status,
values_from = c(
Mean_SC,
SE_SC
)
) %>%
mutate(
Drought_Penalty =
Mean_SC_Drought - Mean_SC_Watered,
Penalty_SE = sqrt(
SE_SC_Drought^2 +
SE_SC_Watered^2
)
)
drought_penalty_data
## # A tibble: 12 × 9
## Date Architecture Mature_Position Mean_SC_Watered Mean_SC_Drought
## <date> <fct> <fct> <dbl> <dbl>
## 1 2026-07-29 Erect Lowest Mature 145. 115.
## 2 2026-07-29 Erect Highest Mature 255. 147.
## 3 2026-07-29 Compact Lowest Mature 200. 130.
## 4 2026-07-29 Compact Highest Mature 206. 178.
## 5 2026-07-30 Erect Lowest Mature 80.8 129.
## 6 2026-07-30 Erect Highest Mature 275. 383.
## 7 2026-07-30 Compact Lowest Mature 122. 158.
## 8 2026-07-30 Compact Highest Mature 241. 275.
## 9 2026-07-31 Erect Lowest Mature 60.3 75.0
## 10 2026-07-31 Erect Highest Mature 162. 162.
## 11 2026-07-31 Compact Lowest Mature 75.6 97.4
## 12 2026-07-31 Compact Highest Mature 159. 144.
## # ℹ 4 more variables: SE_SC_Watered <dbl>, SE_SC_Drought <dbl>,
## # Drought_Penalty <dbl>, Penalty_SE <dbl>
drought_penalty_plot <- ggplot(
drought_penalty_data,
aes(
x = Date,
y = Drought_Penalty,
color = Architecture,
linetype = Mature_Position,
shape = Mature_Position,
group = interaction(
Architecture,
Mature_Position
)
)
) +
geom_hline(
yintercept = 0,
linetype = "dashed",
linewidth = 0.7
) +
geom_errorbar(
aes(
ymin = Drought_Penalty - Penalty_SE,
ymax = Drought_Penalty + Penalty_SE
),
width = 0.08,
linewidth = 0.6,
alpha = 0.75
) +
geom_line(
linewidth = 1.1
) +
geom_point(
size = 3.2
) +
scale_x_date(
date_labels = "%m/%d",
date_breaks = "1 day"
) +
labs(
title = "Architecture Effects on Drought Penalty Across Shoot Positions",
subtitle = "Drought penalty = drought mean − corresponding watered mean",
x = "Date",
y = expression(
Delta *
" Mean Stomatal Conductance (Drought - Watered)"
),
color = "Architecture",
linetype = "Mature Position",
shape = "Mature Position"
) +
theme_classic(
base_size = 12
) +
theme(
plot.title = element_text(
size = 15,
face = "bold",
hjust = 0.5
),
plot.subtitle = element_text(
size = 10,
hjust = 0.5
),
axis.title = element_text(
size = 11,
face = "bold"
),
axis.text = element_text(
size = 10
),
legend.position = "bottom"
)
drought_penalty_plot
##Stability
stability_data <- analysis_data %>%
group_by(
Plant_ID,
Architecture,
Water_Status,
Mature_Position
) %>%
summarise(
n_dates = sum(!is.na(SC)),
Mean_SC = mean(
SC,
na.rm = TRUE
),
SD_SC = sd(
SC,
na.rm = TRUE
),
CV_SC = SD_SC / Mean_SC,
SC_Range =
max(SC, na.rm = TRUE) -
min(SC, na.rm = TRUE),
.groups = "drop"
) %>%
filter(
n_dates >= 2,
is.finite(CV_SC)
)
stability_summary <- stability_data %>%
group_by(
Architecture,
Water_Status,
Mature_Position
) %>%
summarise(
n = n(),
Mean_CV = mean(
CV_SC,
na.rm = TRUE
),
SD_CV = sd(
CV_SC,
na.rm = TRUE
),
SE_CV = SD_CV / sqrt(n),
.groups = "drop"
)
stability_summary
## # A tibble: 8 × 7
## Architecture Water_Status Mature_Position n Mean_CV SD_CV SE_CV
## <fct> <fct> <fct> <int> <dbl> <dbl> <dbl>
## 1 Erect Watered Lowest Mature 4 0.431 0.203 0.102
## 2 Erect Watered Highest Mature 4 0.374 0.147 0.0735
## 3 Erect Drought Lowest Mature 3 0.323 0.128 0.0738
## 4 Erect Drought Highest Mature 3 0.604 0.296 0.171
## 5 Compact Watered Lowest Mature 3 0.506 0.186 0.107
## 6 Compact Watered Highest Mature 3 0.352 0.208 0.120
## 7 Compact Drought Lowest Mature 4 0.273 0.126 0.0632
## 8 Compact Drought Highest Mature 4 0.382 0.283 0.142
stability_plot <- ggplot(
stability_data,
aes(
x = Architecture,
y = CV_SC,
fill = Water_Status
)
) +
geom_boxplot(
position = position_dodge(
width = 0.75
),
width = 0.62,
alpha = 0.65,
outlier.shape = NA
) +
geom_point(
aes(
shape = Mature_Position
),
position = position_jitterdodge(
jitter.width = 0.10,
dodge.width = 0.75
),
size = 2.8,
alpha = 0.85
) +
labs(
title = "Temporal Stability of Stomatal Conductance by Architecture",
subtitle = "Lower coefficients of variation indicate greater stability across dates",
x = "Architecture",
y = "Coefficient of Variation Across Dates",
fill = "Water Status",
shape = "Mature Position"
) +
theme_classic(
base_size = 12
) +
theme(
plot.title = element_text(
size = 15,
face = "bold",
hjust = 0.5
),
plot.subtitle = element_text(
size = 10,
hjust = 0.5
),
axis.title = element_text(
size = 11,
face = "bold"
),
axis.text = element_text(
size = 10
),
legend.position = "bottom"
)
stability_plot
### Interpretation Lower variability indicates more stable conductance
over time. This should be interpreted with the trajectory graph because
a consistently low conductance value can also produce low
variability.
##Architecture × mature-position interaction
interaction_summary <- analysis_data %>%
group_by(
Date,
Architecture,
Water_Status,
Mature_Position
) %>%
summarise(
n = sum(!is.na(SC)),
Mean_SC = mean(
SC,
na.rm = TRUE
),
SD_SC = sd(
SC,
na.rm = TRUE
),
SE_SC = SD_SC / sqrt(n),
.groups = "drop"
)
architecture_position_plot <- ggplot(
interaction_summary,
aes(
x = Mature_Position,
y = Mean_SC,
color = Architecture,
linetype = Water_Status,
shape = Water_Status,
group = interaction(
Architecture,
Water_Status
)
)
) +
geom_line(
linewidth = 1.1,
position = position_dodge(
width = 0.10
)
) +
geom_point(
size = 3.2,
position = position_dodge(
width = 0.10
)
) +
geom_errorbar(
aes(
ymin = Mean_SC - SE_SC,
ymax = Mean_SC + SE_SC
),
width = 0.08,
linewidth = 0.7,
position = position_dodge(
width = 0.10
)
) +
facet_wrap(
~Date,
nrow = 1,
labeller = labeller(
Date = function(x) {
format(
as.Date(x),
"%m/%d"
)
}
)
) +
labs(
title = "Interaction Between Shoot Architecture and Mature Position",
subtitle = "Nonparallel lines suggest architecture affects shoot positions differently",
x = "Mature Meristem Position",
y = expression(
"Mean Stomatal Conductance (mmol " *
m^-2 * " " * s^-1 * ")"
),
color = "Architecture",
linetype = "Water Status",
shape = "Water Status"
) +
theme_classic(
base_size = 12
) +
theme(
plot.title = element_text(
size = 15,
face = "bold",
hjust = 0.5
),
plot.subtitle = element_text(
size = 10,
hjust = 0.5
),
strip.text = element_text(
size = 11,
face = "bold"
),
axis.title = element_text(
size = 11,
face = "bold"
),
axis.text = element_text(
size = 10
),
legend.position = "bottom"
)
architecture_position_plot
redistribution_data <- data_clean %>%
mutate(
Delta_SC =
Highest_Mature_Meristem_SC -
Lowest_Mature_Meristem_SC
)
redistribution_summary <- redistribution_data %>%
group_by(
Date,
Architecture,
Water_Status
) %>%
summarise(
n = n(),
Mean_Delta = mean(
Delta_SC,
na.rm = TRUE
),
SD_Delta = sd(
Delta_SC,
na.rm = TRUE
),
SE_Delta =
SD_Delta / sqrt(n),
.groups = "drop"
)
redistribution_summary
## # A tibble: 12 × 7
## Date Architecture Water_Status n Mean_Delta SD_Delta SE_Delta
## <date> <fct> <fct> <int> <dbl> <dbl> <dbl>
## 1 2026-07-29 Erect Watered 4 110. 31.1 15.6
## 2 2026-07-29 Erect Drought 3 32.2 25.2 14.5
## 3 2026-07-29 Compact Watered 3 6.13 37.5 21.6
## 4 2026-07-29 Compact Drought 4 48.2 60.4 30.2
## 5 2026-07-30 Erect Watered 4 194. 84.7 42.3
## 6 2026-07-30 Erect Drought 3 254. 107. 61.6
## 7 2026-07-30 Compact Watered 3 120. 47.7 27.5
## 8 2026-07-30 Compact Drought 4 117. 87.8 43.9
## 9 2026-07-31 Erect Watered 4 102. 28.1 14.1
## 10 2026-07-31 Erect Drought 3 87.2 63.7 36.8
## 11 2026-07-31 Compact Watered 3 83.5 13.6 7.84
## 12 2026-07-31 Compact Drought 4 46.4 35.3 17.6
redistribution_plot <- ggplot(
redistribution_summary,
aes(
x = Date,
y = Mean_Delta,
color = Architecture,
group = Architecture
)
) +
geom_hline(
yintercept = 0,
linetype = "dashed",
linewidth = 0.7
) +
geom_line(
linewidth = 1.2
) +
geom_point(
size = 3
) +
geom_errorbar(
aes(
ymin = Mean_Delta - SE_Delta,
ymax = Mean_Delta + SE_Delta
),
width = 0.08
) +
facet_wrap(
~Water_Status
) +
scale_x_date(
date_labels = "%m/%d",
date_breaks = "1 day"
) +
labs(
title = "Developmental Redistribution of Stomatal Conductance",
subtitle = "Positive values indicate greater conductance in the highest mature meristem",
x = "Date",
y = expression(
Delta*" Stomatal Conductance (Highest - Lowest)"
),
color = "Architecture"
) +
theme_classic(
base_size = 12
) +
theme(
plot.title = element_text(
face = "bold",
hjust = 0.5
),
plot.subtitle = element_text(
hjust = 0.5
),
strip.text = element_text(
face = "bold"
),
legend.position = "bottom"
)
redistribution_plot
install.packages("lme4")
## Installing package into '/cloud/lib/x86_64-pc-linux-gnu-library/4.6'
## (as 'lib' is unspecified)
install.packages("lmerTest")
## Installing package into '/cloud/lib/x86_64-pc-linux-gnu-library/4.6'
## (as 'lib' is unspecified)
library(lme4)
## Loading required package: Matrix
##
## Attaching package: 'Matrix'
## The following objects are masked from 'package:tidyr':
##
## expand, pack, unpack
library(lmerTest)
##
## Attaching package: 'lmerTest'
## The following object is masked from 'package:lme4':
##
## lmer
## The following object is masked from 'package:stats':
##
## step
redistribution_model <- lmer(
Delta_SC ~
Architecture *
Water_Status *
Date +
(1 | Plant_ID),
data = redistribution_data
)
## Warning: Some predictor variables are on very different scales: consider rescaling.
## You may also use (g)lmerControl(autoscale = TRUE) to improve numerical stability.
## boundary (singular) fit: see help('isSingular')
## Warning: Some predictor variables are on very different scales: consider rescaling.
## You may also use (g)lmerControl(autoscale = TRUE) to improve numerical stability.
anova(redistribution_model)
## Type III Analysis of Variance Table with Satterthwaite's method
## Sum Sq Mean Sq NumDF DenDF F value Pr(>F)
## Date 9974.2 9974.2 1 33.955 1.5077 0.2279
## Architecture:Date 4996.8 4996.8 1 34.056 0.7553 0.3909
## Water_Status:Date 34.6 34.6 1 33.938 0.0052 0.9428
## Architecture:Water_Status:Date 8926.4 8926.4 1 34.013 1.3493 0.2535
## Architecture
## Water_Status
## Architecture:Water_Status