##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