################################################
# Path models - Excluded Traits
################################################
model_data <- read_excel("Output/Data_for_modelling/model_data_v03.xlsx")
model_data$CoverShrubLayer_scaled <- as.vector(scale(asin(sqrt(model_data$CoverShrubLayer))))
model_data$TSF_scaled <- as.vector(scale(asin(sqrt(model_data$TSF))))
model_data$ses_scaled <- as.vector(scale(model_data$ses))
model_data$canopy_num <- as.numeric(ifelse(model_data$canopy == "shady", 0, 1))
model_data$fence_num <- as.numeric(ifelse(model_data$fence == "control", 1, 0))
model_data$pH_scaled <- as.numeric(scale(model_data$pH))
model_data$ses_Leaf_Area_scaled <- as.vector(scale(model_data$ses_Leaf_Area))
model_data$ses_SLA_scaled <- as.vector(scale(model_data$ses_SLA))
model_data$ses_LDMC_scaled <- as.vector(scale(model_data$ses_LDMC))
model_data$ses_Nitrogen_mass_scaled <- as.vector(scale(model_data$ses_Nitrogen_mass))
model_data$ses_Seed_Dry_Mass_scaled <- as.vector(scale(model_data$ses_Seed_Dry_Mass))
model_data$ses_Fruit_Type_scaled <- as.vector(scale(model_data$ses_Fruit_Type))
model_data$ses_Spinescence_scaled <- as.vector(scale(model_data$ses_Spinescence))
model_data$ses_Raunkiaer_LifeForm_scaled <- as.vector(scale(model_data$ses_Raunkiaer_LifeForm))
model_data$ses_Life_Span_scaled <- as.vector(scale(model_data$ses_Life_Span))
model_data$ses_Functional_Type_scaled <- as.vector(scale(model_data$ses_Functional_Type))
model_data$ses_Height_Classification_scaled <- as.vector(scale(model_data$ses_Height_Classification))
model_data$ses_Light_Ellenberg_scaled <- as.vector(scale(model_data$ses_Light_Ellenberg))
model_data$ses_Moisture_Ellenberg_scaled <- as.vector(scale(model_data$ses_Moisture_Ellenberg))
model_data$ses_Nitrogen_Ellenberg_scaled <- as.vector(scale(model_data$ses_Nitrogen_Ellenberg))
model_data$ses_Temperature_Ellenberg_scaled <- as.vector(scale(model_data$ses_Temperature_Ellenberg))
model_data$ses_C_Score_scaled <- as.vector(scale(model_data$ses_C_Score))
model_data$ses_S_Score_scaled <- as.vector(scale(model_data$ses_S_Score))
model_data$ses_R_Score_scaled <- as.vector(scale(model_data$ses_R_Score))
str(model_data)
names(model_data)

sum(is.na(model_data)) # 1 plot has to be removed
model_data <- na.omit(model_data)
sum(is.na(model_data))

# SEM ----
### Shrubs ----
sem.shrub<- lmer(CoverShrubLayer_scaled ~
                   canopy_num +
                   fence_num +
                   pH_scaled +
                   (1|plot),
                 data = model_data)

### TSF ----
sem.TSF<- lmer(TSF_scaled ~ 
                 canopy_num +
                 CoverShrubLayer_scaled + 
                 (1|plot),
               data = model_data)

### ses_Leaf_Area_scaled ----
sem.SESMPD_Leaf_Area <- lmer(ses_Leaf_Area_scaled ~ 
                               TSF_scaled +
                               fence_num +
                               CoverShrubLayer_scaled +
                               pH_scaled +
                               (1|plot), 
                             data = model_data)

## Combining models
model_Leaf_Area = psem(
  sem.shrub,
  sem.TSF,
  sem.SESMPD_Leaf_Area)
x <- summary(model_Leaf_Area)
# Chi-Squared = 2.795 with P-value = 0.424 and on 3 degrees of freedom
# Fisher's C = 13.696 with P-value = 0.033 and on 6 degrees of freedom
x$coefficients[c(6:9),]
Leaf_Area_Light <- x$coefficients[6,8] # signifcant
Leaf_Area_Fence <- x$coefficients[7,8] # 
Leaf_Area_Shrub <- x$coefficients[8,8]
Leaf_Area_pH <- x$coefficients[9,8]

### ses_SLA_scaled ----
sem.SESMPD_SLA <- lmer(ses_SLA_scaled ~ 
                         TSF_scaled +
                         fence_num +
                         CoverShrubLayer_scaled +
                         pH_scaled +
                         (1|plot), 
                       data = model_data)

## Combining models
model_SLA = psem(
  sem.shrub,
  sem.TSF,
  sem.SESMPD_SLA)
x <- summary(model_SLA)
# Chi-Squared = 2.575 with P-value = 0.462 and on 3 degrees of freedom
# Fisher's C = 12.941 with P-value = 0.044 and on 6 degrees of freedom

x$coefficients[c(6:9),]
SLA_Light <- x$coefficients[6,8]
SLA_Fence <- x$coefficients[7,8]
SLA_Shrub <- x$coefficients[8,8]
SLA_pH <- x$coefficients[9,8]

### ses_LDMC_scaled ----
sem.SESMPD_LDMC <- lmer(ses_LDMC_scaled ~ 
                          TSF_scaled +
                          fence_num +
                          CoverShrubLayer_scaled +
                          pH_scaled +
                          (1|plot), 
                        data = model_data)

## Combining models
model_LDMC = psem(
  sem.shrub,
  sem.TSF,
  sem.SESMPD_LDMC)
x <- summary(model_LDMC)
# Chi-Squared = 2.575 with P-value = 0.462 and on 3 degrees of freedom
# Fisher's C = 12.941 with P-value = 0.044 and on 6 degrees of freedom
x$coefficients[c(6:9),]
LDMC_Light <- x$coefficients[6,8] # 
LDMC_Fence <- x$coefficients[7,8] # 
LDMC_Shrub <- x$coefficients[8,8]
LDMC_pH <- x$coefficients[9,8]

### ses_Nitrogen_mass_scaled ----
sem.SESMPD_Nitrogen_mass <- lmer(ses_Nitrogen_mass_scaled ~ 
                                   TSF_scaled +
                                   fence_num +
                                   CoverShrubLayer_scaled +
                                   pH_scaled +
                                   (1|plot), 
                                 data = model_data)

## Combining models
model_Nitrogen_mass = psem(
  sem.shrub,
  sem.TSF,
  sem.SESMPD_Nitrogen_mass)
x <- summary(model_Nitrogen_mass)
# Chi-Squared = 2.577 with P-value = 0.462 and on 3 degrees of freedom
# Fisher's C = 12.897 with P-value = 0.045 and on 6 degrees of freedom
x$coefficients[c(6:9),]
Nitrogen_mass_Light <- x$coefficients[6,8] # 
Nitrogen_mass_Fence <- x$coefficients[7,8] # 
Nitrogen_mass_Shrub <- x$coefficients[8,8]
Nitrogen_mass_pH <- x$coefficients[9,8]

### ses_Seed_Dry_Mass_scaled ----
sem.SESMPD_Seed_Dry_Mass <- lmer(ses_Seed_Dry_Mass_scaled ~ 
                                   TSF_scaled +
                                   fence_num +
                                   CoverShrubLayer_scaled +
                                   pH_scaled +
                                   (1|plot), 
                                 data = model_data)

## Combining models
model_Seed_Dry_Mass = psem(
  sem.shrub,
  sem.TSF,
  sem.SESMPD_Seed_Dry_Mass)
x <- summary(model_Seed_Dry_Mass)
# Chi-Squared = 2.716 with P-value = 0.438 and on 3 degrees of freedom
# Fisher's C = 13.508 with P-value = 0.036 and on 6 degrees of freedom
x$coefficients[c(6:9),]
Seed_Dry_Mass_Light <- x$coefficients[6,8] # 
Seed_Dry_Mass_Fence <- x$coefficients[7,8] # 
Seed_Dry_Mass_Shrub <- x$coefficients[8,8]
Seed_Dry_Mass_pH <- x$coefficients[9,8]

### ses_Fruit_Type_scaled ----
sem.SESMPD_Fruit_Type <- lmer(ses_Fruit_Type_scaled ~ 
                                TSF_scaled +
                                fence_num +
                                CoverShrubLayer_scaled +
                                pH_scaled +
                                (1|plot), 
                              data = model_data)

## Combining models
model_Fruit_Type = psem(
  sem.shrub,
  sem.TSF,
  sem.SESMPD_Fruit_Type)
x <- summary(model_Fruit_Type)
# Chi-Squared = 3.586 with P-value = 0.31 and on 3 degrees of freedom
# Fisher's C = 15.065 with P-value = 0.02 and on 6 degrees of freedom
x$coefficients[c(6:9),]
Fruit_Type_Light <- x$coefficients[6,8] # 
Fruit_Type_Fence <- x$coefficients[7,8] #  
Fruit_Type_Shrub <- x$coefficients[8,8]
Fruit_Type_pH <- x$coefficients[9,8]


### ses_Spinescence_scaled ----
sem.SESMPD_Spinescence <- lmer(ses_Spinescence_scaled ~ 
                                 TSF_scaled +
                                 fence_num +
                                 CoverShrubLayer_scaled +
                                 pH_scaled +
                                 (1|plot), 
                               data = model_data)

## Combining models
model_Spinescence = psem(
  sem.shrub,
  sem.TSF,
  sem.SESMPD_Spinescence)
x <- summary(model_Spinescence)
# Chi-Squared = 2.882 with P-value = 0.41 and on 3 degrees of freedom
# Fisher's C = 13.834 with P-value = 0.032 and on 6 degrees of freedom
x$coefficients[c(6:9),]
Spinescence_Light <- x$coefficients[6,8] # 
Spinescence_Fence <- x$coefficients[7,8] # 
Spinescence_Shrub <- x$coefficients[8,8]
Spinescence_pH <- x$coefficients[9,8]

### ses_Raunkiaer_LifeForm_scaled ----
sem.SESMPD_Raunkiaer_LifeForm <- lmer(ses_Raunkiaer_LifeForm_scaled ~ 
                                        TSF_scaled +
                                        fence_num +
                                        CoverShrubLayer_scaled +
                                        pH_scaled +
                                        (1|plot), 
                                      data = model_data)

## Combining models
model_Raunkiaer_LifeForm = psem(
  sem.shrub,
  sem.TSF,
  sem.SESMPD_Raunkiaer_LifeForm)
x <- summary(model_Raunkiaer_LifeForm)
# Chi-Squared = 2.576 with P-value = 0.462 and on 3 degrees of freedom
# Fisher's C = 12.777 with P-value = 0.047 and on 6 degrees of freedom
x$coefficients[c(6:9),]
Raunkiaer_LifeForm_Light <- x$coefficients[6,8] # 
Raunkiaer_LifeForm_Fence <- x$coefficients[7,8] # 
Raunkiaer_LifeForm_Shrub <- x$coefficients[8,8]
Raunkiaer_LifeForm_pH <- x$coefficients[9,8]

### ses_Life_Span_scaled ----
sem.SESMPD_Life_Span <- lmer(ses_Life_Span_scaled ~ 
                               TSF_scaled +
                               fence_num +
                               CoverShrubLayer_scaled +
                               pH_scaled +
                               (1|plot), 
                             data = model_data)

## Combining models
model_Life_Span = psem(
  sem.shrub,
  sem.TSF,
  sem.SESMPD_Life_Span)
x <- summary(model_Life_Span)
# Chi-Squared = 2.589 with P-value = 0.459 and on 3 degrees of freedom
# Fisher's C = 13.089 with P-value = 0.042 and on 6 degrees of freedom
x$coefficients[c(6:9),]
Life_Span_Light <- x$coefficients[6,8] # 
Life_Span_Fence <- x$coefficients[7,8] # 
Life_Span_Shrub <- x$coefficients[8,8]
Life_Span_pH <- x$coefficients[9,8]


### ses_Functional_Type_scaled ----
sem.SESMPD_Functional_Type <- lmer(ses_Functional_Type_scaled ~ 
                                     TSF_scaled +
                                     fence_num +
                                     CoverShrubLayer_scaled +
                                     pH_scaled +
                                     (1|plot), 
                                   data = model_data)

## Combining models
model_Functional_Type = psem(
  sem.shrub,
  sem.TSF,
  sem.SESMPD_Functional_Type)
x <- summary(model_Functional_Type)
# Chi-Squared = 2.655 with P-value = 0.448 and on 3 degrees of freedom
# Fisher's C = 13.331 with P-value = 0.038 and on 6 degrees of freedom
x$coefficients[c(6:9),]
Functional_Type_Light <- x$coefficients[6,8] # 
Functional_Type_Fence <- x$coefficients[7,8] # 
Functional_Type_Shrub <- x$coefficients[8,8]
Functional_Type_pH <- x$coefficients[9,8]

### ses_Height_Classification_scaled ----
sem.SESMPD_Height_Classification <- lmer(ses_Height_Classification_scaled ~ 
                                           TSF_scaled +
                                           fence_num +
                                           CoverShrubLayer_scaled +
                                           pH_scaled +
                                           (1|plot), 
                                         data = model_data)

## Combining models
model_Height_Classification = psem(
  sem.shrub,
  sem.TSF,
  sem.SESMPD_Height_Classification)
x <- summary(model_Height_Classification)
# Chi-Squared = 2.627 with P-value = 0.453 and on 3 degrees of freedom
# Fisher's C = 13.246 with P-value = 0.039 and on 6 degrees of freedom
x$coefficients[c(6:9),]
Height_Classification_Light <- x$coefficients[6,8] # 
Height_Classification_Fence <- x$coefficients[7,8] # 
Height_Classification_Shrub <- x$coefficients[8,8]
Height_Classification_pH <- x$coefficients[9,8]

### ses_Light_Ellenberg_scaled ----
sem.SESMPD_Light_Ellenberg <- lmer(ses_Light_Ellenberg_scaled ~ 
                                     TSF_scaled +
                                     fence_num +
                                     CoverShrubLayer_scaled +
                                     pH_scaled +
                                     (1|plot), 
                                   data = model_data)

## Combining models
model_Light_Ellenberg = psem(
  sem.shrub,
  sem.TSF,
  sem.SESMPD_Light_Ellenberg)
x <- summary(model_Light_Ellenberg)
# Chi-Squared = 2.66 with P-value = 0.447 and on 3 degrees of freedom
# Fisher's C = 13.304 with P-value = 0.038 and on 6 degrees of freedom
x$coefficients[c(6:9),]
Light_Ellenberg_Light <- x$coefficients[6,8] # 
Light_Ellenberg_Fence <- x$coefficients[7,8] # 
Light_Ellenberg_Shrub <- x$coefficients[8,8]
Light_Ellenberg_pH <- x$coefficients[9,8]

### ses_Moisture_Ellenberg_scaled ----
sem.SESMPD_Moisture_Ellenberg <- lmer(ses_Moisture_Ellenberg_scaled ~ 
                                        TSF_scaled +
                                        fence_num +
                                        CoverShrubLayer_scaled +
                                        pH_scaled +
                                        (1|plot), 
                                      data = model_data)

## Combining models
model_Moisture_Ellenberg = psem(
  sem.shrub,
  sem.TSF,
  sem.SESMPD_Moisture_Ellenberg)
x <- summary(model_Moisture_Ellenberg)
# Chi-Squared = 2.641 with P-value = 0.45 and on 3 degrees of freedom
# Fisher's C = 13.274 with P-value = 0.039 and on 6 degrees of freedom
x$coefficients[c(6:9),]
Moisture_Ellenberg_Light <- x$coefficients[6,8] # 
Moisture_Ellenberg_Fence <- x$coefficients[7,8] # signifcant
Moisture_Ellenberg_Shrub <- x$coefficients[8,8]
Moisture_Ellenberg_pH <- x$coefficients[9,8]

### ses_Nitrogen_Ellenberg_scaled ----
sem.SESMPD_Nitrogen_Ellenberg <- lmer(ses_Nitrogen_Ellenberg_scaled ~ 
                                        TSF_scaled +
                                        fence_num +
                                        CoverShrubLayer_scaled +
                                        pH_scaled +
                                        (1|plot), 
                                      data = model_data)

## Combining models
model_Nitrogen_Ellenberg = psem(
  sem.shrub,
  sem.TSF,
  sem.SESMPD_Nitrogen_Ellenberg)
x <- summary(model_Nitrogen_Ellenberg)
# Chi-Squared = 2.991 with P-value = 0.393 and on 3 degrees of freedom
# Fisher's C = 14.132 with P-value = 0.028 and on 6 degrees of freedom
x$coefficients[c(6:9),]
Nitrogen_Ellenberg_Light <- x$coefficients[6,8] # 
Nitrogen_Ellenberg_Fence <- x$coefficients[7,8] # 
Nitrogen_Ellenberg_Shrub <- x$coefficients[8,8]
Nitrogen_Ellenberg_pH <- x$coefficients[9,8]

### ses_Temperature_Ellenberg_scaled ----
sem.SESMPD_Temperature_Ellenberg <- lmer(ses_Temperature_Ellenberg_scaled ~ 
                                           TSF_scaled +
                                           fence_num +
                                           CoverShrubLayer_scaled +
                                           pH_scaled +
                                           (1|plot), 
                                         data = model_data)

## Combining models
model_Temperature_Ellenberg = psem(
  sem.shrub,
  sem.TSF,
  sem.SESMPD_Temperature_Ellenberg)
x <- summary(model_Temperature_Ellenberg)
# Chi-Squared = 3.096 with P-value = 0.377 and on 3 degrees of freedom
# Fisher's C = 14.304 with P-value = 0.026 and on 6 degrees of freedom
x$coefficients[c(6:9),]
Temperature_Ellenberg_Light <- x$coefficients[6,8] # 
Temperature_Ellenberg_Fence <- x$coefficients[7,8] # 
Temperature_Ellenberg_Shrub <- x$coefficients[8,8]
Temperature_Ellenberg_pH <- x$coefficients[9,8]

### ses_C_Score_scaled ----
sem.SESMPD_C_Score <- lmer(ses_C_Score_scaled ~ 
                             TSF_scaled +
                             fence_num +
                             CoverShrubLayer_scaled +
                             pH_scaled +
                             (1|plot), 
                           data = model_data)

## Combining models
model_C_Score = psem(
  sem.shrub,
  sem.TSF,
  sem.SESMPD_C_Score)
x <- summary(model_C_Score)
# Chi-Squared = 2.615 with P-value = 0.455 and on 3 degrees of freedom
# Fisher's C = 13.149 with P-value = 0.041 and on 6 degrees of freedom
x$coefficients[c(6:9),]
C_Score_Light <- x$coefficients[6,8] # 
C_Score_Fence <- x$coefficients[7,8] # 
C_Score_Shrub <- x$coefficients[8,8]
C_Score_pH <- x$coefficients[9,8]

### ses_S_Score_scaled ----
sem.SESMPD_S_Score <- lmer(ses_S_Score_scaled ~ 
                             TSF_scaled +
                             fence_num +
                             CoverShrubLayer_scaled +
                             pH_scaled +
                             (1|plot), 
                           data = model_data)

## Combining models
model_S_Score = psem(
  sem.shrub,
  sem.TSF,
  sem.SESMPD_S_Score)
x <- summary(model_S_Score)
# Chi-Squared = 3.006 with P-value = 0.391 and on 3 degrees of freedom
# Fisher's C = 14.175 with P-value = 0.028 and on 6 degrees of freedom
x$coefficients[c(6:9),]
S_Score_Light <- x$coefficients[6,8] # 
S_Score_Fence <- x$coefficients[7,8] # 
S_Score_Shrub <- x$coefficients[8,8]
S_Score_pH <- x$coefficients[9,8]

### ses_R_Score_scaled ----
sem.SESMPD_R_Score <- lmer(ses_R_Score_scaled ~ 
                             TSF_scaled +
                             fence_num +
                             CoverShrubLayer_scaled +
                             pH_scaled +
                             (1|plot), 
                           data = model_data)

## Combining models
model_R_Score = psem(
  sem.shrub,
  sem.TSF,
  sem.SESMPD_R_Score)
x <- summary(model_R_Score)
# Chi-Squared = 2.638 with P-value = 0.451 and on 3 degrees of freedom
# Fisher's C = 13.302 with P-value = 0.038 and on 6 degrees of freedom
x$coefficients[c(6:9),]
R_Score_Light <- x$coefficients[6,8] # 
R_Score_Fence <- x$coefficients[7,8] # 
R_Score_Shrub <- x$coefficients[8,8]
R_Score_pH <- x$coefficients[9,8]



## summarize ----
# Create a data frame summarizing all extracted coefficients
coef_table <- data.frame(
  Trait = c(
    "Leaf area", "Specific leaf area", "Leaf dry matter content", "Leaf nitrogen content",
    "Seed mass", "Fruit type", "Spinescence",
    "Raunkiaer plant life form", "Plant lifespan", "Plant functional type",
    "Plant height", "Ellenberg light", 
    "Ellenberg moisture", "Ellenberg nitrogen", 
    "Ellenberg temperature", "Grime's C-score", "Grime's S-score", "Grime's R-score"
  ),
  Light = c(
    Leaf_Area_Light, SLA_Light, LDMC_Light, Nitrogen_mass_Light,
    Seed_Dry_Mass_Light, Fruit_Type_Light, Spinescence_Light,
    Raunkiaer_LifeForm_Light, Life_Span_Light, Functional_Type_Light,
    Height_Classification_Light, Light_Ellenberg_Light,
    Moisture_Ellenberg_Light, Nitrogen_Ellenberg_Light,
    Temperature_Ellenberg_Light, C_Score_Light, S_Score_Light, R_Score_Light
  ),
  Fence = c(
    Leaf_Area_Fence, SLA_Fence, LDMC_Fence, Nitrogen_mass_Fence,
    Seed_Dry_Mass_Fence, Fruit_Type_Fence, Spinescence_Fence,
    Raunkiaer_LifeForm_Fence, Life_Span_Fence, Functional_Type_Fence,
    Height_Classification_Fence, Light_Ellenberg_Fence,
    Moisture_Ellenberg_Fence, Nitrogen_Ellenberg_Fence,
    Temperature_Ellenberg_Fence, C_Score_Fence, S_Score_Fence, R_Score_Fence
  ),
  Shrub = c(
    Leaf_Area_Shrub, SLA_Shrub, LDMC_Shrub, Nitrogen_mass_Shrub,
    Seed_Dry_Mass_Shrub, Fruit_Type_Shrub, Spinescence_Shrub,
    Raunkiaer_LifeForm_Shrub, Life_Span_Shrub, Functional_Type_Shrub,
    Height_Classification_Shrub, Light_Ellenberg_Shrub,
    Moisture_Ellenberg_Shrub, Nitrogen_Ellenberg_Shrub,
    Temperature_Ellenberg_Shrub, C_Score_Shrub, S_Score_Shrub, S_Score_Shrub
  ),
  pH = c(
    Leaf_Area_pH, SLA_pH, LDMC_pH, Nitrogen_mass_pH,
    Seed_Dry_Mass_pH, Fruit_Type_pH, Spinescence_pH,
    Raunkiaer_LifeForm_pH, Life_Span_pH, Functional_Type_pH,
    Height_Classification_pH, Light_Ellenberg_pH,
    Moisture_Ellenberg_pH, Nitrogen_Ellenberg_pH,
    Temperature_Ellenberg_pH, C_Score_pH, S_Score_pH, S_Score_Shrub
  )
)

graphics.off()

png("Output/Grafics/Excluded_Traits.png", 
    width = 220, height = 200, 
    units = "mm", res = 900)
# Ordered data
df_f <- coef_table[order(coef_table$Fence), ]
df_l <- coef_table[order(coef_table$Light), ]
# df_l <- coef_table[order(coef_table$Fence), ]

n <- nrow(df_f)
ypos <- seq(n, 1)   # shared vertical positions
bar_height <- 0.35  # half-height of each bar

# Two panels
par(mfrow = c(1, 2),
    oma   = c(3, 4, 3, 4),
    xaxs  = "i")

##########################################
### LEFT PANEL (Roe deer presence)
##########################################
par(mar = c(2, 5, 0.5, 0))

## Manual x-axis ticks & labels (LEFT)
ticks_left  <- c(-0.25, -0.20, -0.15, -0.10, -0.05, 0)
labels_left <- c("-0.25", "-0.20", "-0.15", "-0.10", "-0.05", "0")

# x-axis limits must include all ticks
xlim_left <- c(min(ticks_left) - 0.01, 0)

plot(NULL,
     xlim = xlim_left,
     ylim = c(0.5, n + 0.5),
     xaxt = "n", yaxt = "n",
     xlab = "", ylab = "")
mtext("Roe deer presence", 
      side = 3, 
      line = 0.5, 
      cex = 1.2, 
      font = 2)

# Draw rectangular bars
for (i in 1:n) {
  rect(xleft  = df_f$Fence[i],
       xright = 0,
       ybottom = ypos[i] - bar_height,
       ytop    = ypos[i] + bar_height,
       col = "#016064",
       border = "black")
}

axis(1, at = ticks_left, labels = labels_left)
axis(2, at = ypos, labels = df_f$Trait, las = 1, cex.axis = 0.8)
box()
abline(v = -0.2050, lwd = 2, lty = 2, col = "black")

##########################################
### RIGHT PANEL (Light availability)
##########################################
par(mar = c(2, 0, 0.5, 5))

## Manual x-axis ticks & labels (RIGHT)
ticks_right  <- c(0, 0.05, 0.10, 0.15, 0.20, 0.25, 0.30)
labels_right <- c("0", "0.05", "0.10", "0.15", "0.20", "0.25", "0.30")

# x-axis limits must include all ticks
xlim_right <- c(0, max(ticks_right) + 0.01)

plot(NULL,
     xlim = xlim_right,
     ylim = c(0.5, n + 0.5),
     xaxt = "n", yaxt = "n",
     xlab = "", ylab = "")
mtext("Light availability", 
      side = 3, 
      line = ,.5, 
      cex = 1.2, 
      font = 2)

# Draw rectangular bars
for (i in 1:n) {
  rect(xleft  = 0,
       xright = df_l$Light[i],
       ybottom = ypos[i] - bar_height,
       ytop    = ypos[i] + bar_height,
       col = "#CC7722",
       border = "black")
}

axis(1, at = ticks_right, labels = labels_right)
axis(4, at = ypos, labels = df_l$Trait, las = 1, cex.axis = 0.8)
box()
abline(v = 0.2443, lwd = 2, lty = 2, col = "black")

##########################################
### Shared x axis label
##########################################
mtext("Standardized estimate", side = 1, outer = TRUE, line = 0.5)

dev.off()

# overall model fence: -0.2050
# overall model light: 0.2443


## full-model standardized paths
full_light <- 0.2443
full_fence <- -0.2050

##
coef_table$Delta_Light <- full_light - coef_table$Light
coef_table$Delta_Fence <- full_fence - coef_table$Fence

coef_table

coef_table_fence_ranked <- coef_table[order(-coef_table$Delta_Fence), ]
coef_table_fence_ranked

coef_table_light_ranked <- coef_table[order(-coef_table$Delta_Light), ]
coef_table_light_ranked