---
title: "CO2_analysis_labexp3"
author: "Arthur Vienne"
date: "2025-12-10"
output: html_document
---

```{r}
folder_path <-  "C:/Users/AVienne/OneDrive - Universiteit Antwerpen/Labexp 3 (start Oct 2022)/Picarro/Slopes"
files <- list.files(path = folder_path, full.names = TRUE)
# Read all the matched Excel files into a list
all_slopefiles <- lapply(files, read.csv)

for (i in (1:32)) {all_slopefiles[[i]]<- all_slopefiles[[i]] [,colnames(all_slopefiles[[1]])] 
} #this df had more columns}

all_slopefiles_df <- do.call(rbind, all_slopefiles)   
all_slopefiles_df$Date <- gsub("19_01_2023", "16_01_2023", all_slopefiles_df$Date)

#we did not measure at 19/1/2023, the file should be 16/1/2023, day 55 is the correct one!!!
#Replace run two sorting scripts again!!!
#Merged in the line above, now sort & add time in days after amendment
library (lubridate)
all_slopefiles_df$day <- as.numeric(dmy(all_slopefiles_df$Date) - dmy("22/11/2022") ) #amendmentdate = subtracted here
 all_slopefiles_df <- all_slopefiles_df [which(all_slopefiles_df$R2_CO2_12>0.9 & all_slopefiles_df$R2_CO2_13>0.9 ),]
 
 #1046 / 1182 , approx. 88% of the data is kept

# plot (all_slopefiles_df$day,all_slopefiles_df$dCH4_ppm_s, ylim = c(0,1))
# points(all_slopefiles_df[which(all_slopefiles_df$R2_CH4<0.9),]$day, all_slopefiles_df[which(all_slopefiles_df$R2_CH4<0.9),]$dCH4_ppm_s, col = "red")

#fill data gaps for temperature from another file to calculate density 
library(readxl)
manual_co2_data <- read_excel("C:/Users/AVienne/OneDrive - Universiteit Antwerpen/Labexp 3 (start Oct 2022)/CO2_CH4_dataframe.xlsx",                        sheet = "CO2dataframe")
indices_NA_temp <- c(which(is.na (all_slopefiles_df$T_.C))  , which (all_slopefiles_df$T_.C==0)) #600 points

for (i in (1: length (indices_NA_temp))){
  index <- indices_NA_temp[i]
  DAY <-   all_slopefiles_df[index, ]$day
  Date <-   all_slopefiles_df[index, ]$Date
  
  if (DAY > max (manual_co2_data$Day)){
    #add the correct temperature for all these dates manually
    if (Date=="24_10_2023") {all_slopefiles_df[index, ]$T_.C <- 18     }
    if (Date=="25_10_2023") {all_slopefiles_df[index, ]$T_.C <- 19     }
    if (Date=="9_11_2023") {all_slopefiles_df[index, ]$T_.C <- 21    }
    if (Date=="23_11_2023") {all_slopefiles_df[index, ]$T_.C <- 20.5    }
    if (Date=="13_12_2023") {all_slopefiles_df[index, ]$T_.C <- 19    }
  }
  else {   
    pot <- all_slopefiles_df[index, ]$PotNo
    all_slopefiles_df[index, ]$T_.C <-     manual_co2_data[which(manual_co2_data$Day == DAY & manual_co2_data$potNo==pot), ]$temperature
  }}



### do the same for moisture
# Indices where SoilMoisture is NA or zero
indices_NA_moisture <- c(
  which(is.na(all_slopefiles_df$SoilMoisture)),
  which(all_slopefiles_df$SoilMoisture == 0)
)

for (i in 1:length(indices_NA_moisture)) {
  index <- indices_NA_moisture[i]
  DAY <- all_slopefiles_df[index, ]$day
  Date <- all_slopefiles_df[index, ]$Date
  pot <- all_slopefiles_df[index, ]$PotNo
  if (DAY == 255) next
  all_slopefiles_df[index, ]$SoilMoisture <- 100*manual_co2_data[which(manual_co2_data$Day == DAY & manual_co2_data$potNo == pot),]$SWC
  }

all_slopefiles_df$temperature <- all_slopefiles_df$T_.C



#do an identical loop for SWC, note that there were some soilmoisture data points as SWC (0.xx) and some as %, filter these out!!
#for statistics of slope and SWC
all_slopefiles_df$CH4density <- 101325*16/(8.314*(273+all_slopefiles_df$T_.C))/1000
all_slopefiles_df$CO2density <- 101325*44/(8.314*(273+all_slopefiles_df$T_.C))/1000

#convert to µmol GHG/m2/s:

#unit check: µLCO2/L air/s *gCH4/L CO2 * L air in chamber / Area chamber / gCO2/mol CO2 = µmol CO2/m2/s, same for CH4
all_slopefiles_df$µmolCH4m2s <- all_slopefiles_df$CH4slope*all_slopefiles_df$CH4density*0.98125/0.008171/16
all_slopefiles_df$µmolCO2m2s <- as.numeric (all_slopefiles_df$CO2slope)*all_slopefiles_df$CO2density*0.98125/0.008171/44

data <- all_slopefiles_df

data$Treatment <- NA  # Initialize the Treatment column
# Assign treatments based on PotNo ranges
data$Treatment[which (data$PotNo >= 1 & data$PotNo <= 5)] <- "C3"
data$Treatment[which(data$PotNo >= 6 & data$PotNo <= 10)] <- "FB3"
data$Treatment[data$PotNo >= 11 & data$PotNo <= 15] <- "CB3"
data$Treatment[data$PotNo >= 16 & data$PotNo <= 20] <- "CD3"
data$Treatment[data$PotNo >= 21 & data$PotNo <= 25] <- "MD3"
data$Treatment[data$PotNo >= 26 & data$PotNo <= 30] <- "FD3"
data$Treatment[data$PotNo >= 31 & data$PotNo <= 35] <- "C1"
data$Treatment[data$PotNo >= 41 & data$PotNo <= 45] <- "FB1"
data$Day <- data$day

#add columns to data 
data$SOC  <- as.numeric (with(data,
                              ifelse(Treatment == "C1", 1,
                                     ifelse(Treatment == "FB1", 1, 3))))

data$Basalt <- as.numeric (with(data,
                                ifelse(Treatment == "FB1", 1,
                                       ifelse(Treatment == "FB3", 1,
                                              ifelse(Treatment == "CB3", 1,0)))))

data$Dunite <- as.numeric (with(data,
                                ifelse(Treatment == "FD3", 1,
                                       ifelse(Treatment == "CD3", 1,
                                              ifelse(Treatment == "MD3", 1,0)))))


data$rock <-  with(data,
                   ifelse(Treatment == "FB3", "Basalt",
                          ifelse(Treatment == "CB3","Basalt" ,
                                 ifelse(Treatment == "FB1","Basalt" ,
                                        ifelse(Treatment == "C1","aNone" ,
                                               ifelse(Treatment == "C3", "aNone","Dunite"))))))


data$sizedunite <- as.numeric (with(data,
                                    ifelse(Treatment == "FD3", 71,
                                           ifelse(Treatment == "MD3", 278,
                                                  ifelse(Treatment == "CD3", 1030,
                                                         ifelse(Treatment %in% c("FB1", "FB3"), 0,
                                                                ifelse(Treatment == "CB3", 0, 0)))))))

data$sizebasalt <- as.numeric (with(data, 
                                    ifelse(Treatment == "FD3", 0, 
                                           ifelse(Treatment == "MD3", 0, 
                                                  ifelse(Treatment == "CD3", 0, 
                                                         ifelse(Treatment %in% c("FB1", "FB3"), 292, 
                                                                ifelse(Treatment == "CB3", 670, 0)))))))
data$treatment <- data$Treatment
data$SWC <- data$SoilMoisture / 100
data$Temperature <- data$temperature

data_original <- data
```


```{r}
#CO2 plot 
dataCO2 <- data[complete.cases(data$µmolCO2m2s),]
# #Normalize for T & SWC
#### plots
dataCO2_1 <- dataCO2[which(dataCO2$SOC==1),]
# Step 1: Summarize data by Day
summary_data <- summarySE(dataCO2_1, 
                          measurevar = "µmolCO2m2s", 
                          groupvars = c("Day", "Treatment"))

#left out visualization because outliers and low replicates
summary_data<- summary_data[-which (summary_data$Day==41 ),]
summary_data<- summary_data[-which (summary_data$Day==23 ),]


# Step 2: Plot the data
CO2_flux_plot_1 <- ggplot(summary_data, aes_string(
  x = "Day", 
  y = "µmolCO2m2s", 
  col = "Treatment", 
  shape = "Treatment", 
  size = "Treatment", 
  fill = "Treatment"
)) +
  theme_bw() +
  theme(
    axis.title.y = element_text(size = 10), 
    axis.title.x = element_text(size = 10), 
    axis.title = element_text(face = "bold"),
    panel.grid.major = element_blank(), 
    panel.grid.minor = element_blank(), 
    panel.background = element_blank(), 
    axis.line = element_line(colour = "black"),
    legend.position = "right"
  ) +
  geom_point() +   geom_line(size = 0.5) +
  geom_errorbar(aes(ymin = µmolCO2m2s - se, ymax = µmolCO2m2s + se), size = 0.1) +
  
  # Match manual aesthetics
  scale_shape_manual(values = c(
    "C1" = 22, "FB1" = 21, "C3" = 22, "CB3" = 21, 
    "CD3" = 23, "FB3" = 21, "FD3" = 23, "MD3" = 23)) +
  
  scale_size_manual(values = c(
    "C1" = 4, "FB1" = 2, "C3" = 4, "CB3" = 6, 
    "CD3" = 6, "FB3" = 2, "FD3" = 2, "MD3" = 4)) +
  
  scale_color_manual(values = c(
    "C1" = "black", "FB1" = "grey", "C3" = "black", "CB3" = "grey", 
    "CD3" = "darkgreen", "FB3" = "grey", "FD3" = "darkgreen", "MD3" = "darkgreen")) +
  
  scale_fill_manual(values = c(
    "C1" = "white", "FB1" = "white", "C3" = "black", "CB3" = "grey", 
    "CD3" = "darkgreen", "FB3" = "grey", "FD3" = "darkgreen", "MD3" = "darkgreen")) +
  
  # Labels
  ylab(bquote("SCE" ~ "\n" ~ (mu * "mol CO"[2] ~ m^{-2} ~ s^{-1})))+
  xlab("Time \n (days)")


#now for 3% soils
dataCO2_3 <- dataCO2[-which(dataCO2$SOC==1),]
summary_data <- summarySE(dataCO2_3, 
                          measurevar = "µmolCO2m2s", 
                          groupvars = c("Day", "Treatment"))

CO2_flux_plot_3 <- ggplot(summary_data, aes_string(
  x = "Day", 
  y = "µmolCO2m2s", 
  col = "Treatment", 
  shape = "Treatment", 
  size = "Treatment", 
  fill = "Treatment"
)) +
  theme_bw() +
  theme(
    axis.title.y = element_text(size = 10), 
    axis.title.x = element_text(size = 10), 
    axis.title = element_text(face = "bold"),
    panel.grid.major = element_blank(), 
    panel.grid.minor = element_blank(), 
    panel.background = element_blank(), 
    axis.line = element_line(colour = "black"),
    legend.position = "right"
  ) +
  geom_point() +   geom_line(size = 0.5) +
  geom_errorbar(aes(ymin = µmolCO2m2s - se, ymax = µmolCO2m2s + se), size = 0.1) +
  
  # Match manual aesthetics
  scale_shape_manual(values = c(
    "C1" = 22, "FB1" = 21, "C3" = 22, "CB3" = 21, 
    "CD3" = 23, "FB3" = 21, "FD3" = 23, "MD3" = 23)) +
  
  scale_size_manual(values = c(
    "C1" = 4, "FB1" = 2, "C3" = 4, "CB3" = 6, 
    "CD3" = 6, "FB3" = 2, "FD3" = 2, "MD3" = 4)) +
  
  scale_color_manual(values = c(
    "C1" = "black", "FB1" = "grey", "C3" = "black", "CB3" = "grey", 
    "CD3" = "darkgreen", "FB3" = "grey", "FD3" = "darkgreen", "MD3" = "darkgreen")) +
  
  scale_fill_manual(values = c(
    "C1" = "white", "FB1" = "white", "C3" = "black", "CB3" = "grey", 
    "CD3" = "darkgreen", "FB3" = "grey", "FD3" = "darkgreen", "MD3" = "darkgreen")) +
  
  # Labels
  ylab(bquote("SCE" ~ "\n" ~ (mu * "mol CO"[2] ~ m^{-2} ~ s^{-1})))+
xlab("Time \n (days)")

# Load cowplot
library(cowplot)

# Combine the plots with labels A and B
combined_plot <- plot_grid(
  CO2_flux_plot_1, 
  CO2_flux_plot_3, 
  labels = c("A", "B"),
  label_size = 14,
  ncol = 1,        # ⬅ vertical alignment
  align = "v"      # align vertically
)


setwd ("C:/Users/AVienne/OneDrive - Universiteit Antwerpen/Labexp 3 (start Oct 2022)/PAPER_")
# Save the combined figure at 600 dpi
ggsave("CO2_flux_combined_plot.jpeg", plot = combined_plot,
       width = 12, height = 12, dpi = 600, units = "in")
```


```{r}
#moisture and temperature plots
T_plot_exp <- ggplot(data, aes(x = day, y = Temperature)) +
  geom_point(color = "steelblue", size = 2) +geom_line()+
  xlab("Day after start of experiment") +
  ylab("Temperature (°C)") +
  theme_minimal(base_size = 14)+
  geom_vline(xintercept = 130, linetype = "dotted", color = "black", size = 1) 

ggsave("temperature.jpeg", plot = T_plot_exp,  dpi = 600, units = "in")



#delete days 247, 255 & 267 from the plot, then no good SWC data
data_original_T_SWC_SWCplot <- data[-which (data$day == 255 |data$day == 247 | data$day == 267),]
forlegend <- ggplot(data_original_T_SWC_SWCplot,
       aes(x = day, y = SWC,
           shape = Treatment, col = Treatment,
           size = Treatment, fill = Treatment)) +
  geom_point(size = 3) +
  xlab("Day after start of experiment") + 
  ylab("Soil moisture (m³/m³)") +
  theme_minimal(base_size = 14) +
  scale_shape_manual(values = c(
    "C1" = 22, "FB1" = 21, "C3" = 22, "CB3" = 21, 
    "CD3" = 23, "FB3" = 21, "FD3" = 23, "MD3" = 23
  )) +
  scale_size_manual(values = c(
    "C1" = 4, "FB1" = 2, "C3" = 4, "CB3" = 6, 
    "CD3" = 6, "FB3" = 2, "FD3" = 2, "MD3" = 4
  )) +
  geom_vline(xintercept = 130, linetype = "dotted",
             color = "black", size = 1) +
  scale_color_manual(values = c(
    "C1" = "black", "FB1" = "grey", "C3" = "black", "CB3" = "grey", 
    "CD3" = "darkgreen", "FB3" = "grey", "FD3" = "darkgreen", "MD3" = "darkgreen"
  )) +
  scale_fill_manual(values = c(
    "C1" = "white", "FB1" = "white", "C3" = "black", "CB3" = "grey", 
    "CD3" = "darkgreen", "FB3" = "grey", "FD3" = "darkgreen", "MD3" = "darkgreen"
  )) +
  theme(
    legend.position = "bottom",
    legend.direction = "horizontal"
  )


SWC_plot_exp<- ggplot(data_original_T_SWC_SWCplot, aes(x = day, y = SWC, shape = Treatment,col = Treatment,size = Treatment, fill = Treatment)) +
  geom_point(size = 3) +  # adjust point size if needed
  xlab("Day after start of experiment") +
  ylab("Soil moisture (m³/m³)") +
  theme_minimal(base_size = 14) +
  scale_shape_manual(values = c(
    "C1" = 22, "FB1" = 21, "C3" = 22, "CB3" = 21, 
    "CD3" = 23, "FB3" = 21, "FD3" = 23, "MD3" = 23
  ))+
   scale_size_manual(values = c(
    "C1" = 4, "FB1" = 2, "C3" = 4, "CB3" = 6, 
    "CD3" = 6, "FB3" = 2, "FD3" = 2, "MD3" = 4)) +
  geom_vline(xintercept = 130, linetype = "dotted", color = "black", size = 1) +
  scale_color_manual(values = c(
    "C1" = "black", "FB1" = "grey", "C3" = "black", "CB3" = "grey", 
    "CD3" = "darkgreen", "FB3" = "grey", "FD3" = "darkgreen", "MD3" = "darkgreen")) +
  scale_fill_manual(values = c(
    "C1" = "white", "FB1" = "white", "C3" = "black", "CB3" = "grey", 
    "CD3" = "darkgreen", "FB3" = "grey", "FD3" = "darkgreen", "MD3" = "darkgreen"))+
  theme(legend.position = "none")
  

ggsave("SWC.jpeg", plot = SWC_plot_exp,  dpi = 600, width = 10, height = 10, units = "in")


library (cowplot)
legend <- get_legend (forlegend)
combined_plot <- plot_grid(
  T_plot_exp, SWC_plot_exp, legend,
  labels = c("A", "B", ""),     # labels for the panels
  ncol = 1, nrow =3, 
  rel_heights = c(2, 2, 1)    # adjust height for the legend
)


# Save the combined figure at 600 dpi
ggsave("temperature_SWC_subplot.jpeg",
       plot = combined_plot,
       width = 12,
       height = 10,
       units = "in",
       dpi = 600)

```

```{r}
#data analysis for CO2: preprocessing
data<- data[complete.cases (data$temperature),]
data<- data[complete.cases (data$SWC),]
data<- data[complete.cases (data$µmolCO2m2s),] 
#only a few datapoints where no CO2 was measured in some pots
data$SWC2 <- data$SWC^2

##divide dataset in two datasets
dataCO2FB_C <- data[which(data$Treatment=="C1"|data$Treatment=="C3"|data$Treatment=="FB3"|data$Treatment=="FB1"),]
dataCO2D_B_size <- data[-which(data$Treatment=="C1"|data$Treatment=="FB1"),]


#for dataCO2FB_C, CO2
#1)normalization
modelCO2 <- lmer (data = dataCO2FB_C, µmolCO2m2s ~ Temperature + SWC + SWC2 +  (1 | PotNo) )
summary (modelCO2)
Anova (modelCO2)
dataCO2FB_C$residualsCO2 <- residuals(modelCO2)

#2 model
modelCO2 <- lmer (data = dataCO2FB_C, residualsCO2 ~ SOC + Basalt + Day + Basalt:Day + Basalt:SOC + (1 | PotNo) )
summary (modelCO2)
Anova (modelCO2)

#N.S. for all interactions, new model without interactions
modelCO2 <- lmer (data = dataCO2FB_C, residualsCO2 ~ SOC + Basalt + Day+ (1 | PotNo) )
summary (modelCO2)
Anova (modelCO2)
qqnorm(residuals(modelCO2))
car::vif(modelCO2) #VIF = 5, still OK. 


#### for dataset 2: 
library(emmeans)
modelCO2 <- lmer (data = dataCO2D_B_size, µmolCO2m2s ~ Temperature + SWC + SWC2  +  (1 | PotNo) )
summary (modelCO2)
Anova (modelCO2)
dataCO2D_B_size$residualsCO2 <- residuals(modelCO2)

#check interactions , rock:day gave p =0.07 N.S.
modelCO2 <- lmer (data = dataCO2D_B_size, residualsCO2 ~   sizebasalt + sizedunite + Day + rock +   (1 | PotNo) )
summary (modelCO2)
Anova (modelCO2)
qqnorm(residuals(modelCO2))
car::vif(modelCO2) #VIF < 5 OK
emmeans(modelCO2, pairwise ~ rock)



```



```{r}
#Q10s:  for all data together
data_original$logCO2 <- log(data_original$µmolCO2m2s)
model_Q10 <- lmer(logCO2 ~ Temperature * rock + (1 | PotNo),
                  data = data_original)
summary(model_Q10)
Anova(model_Q10)
# If Temperature:rock interaction is significant → Q10 differs among rocks
# Extract slopes per rock (b parameter)
slopes <- emtrends(model_Q10, ~ rock, var = "Temperature")
# Pairwise tests: Are Q10 values different?
pairs(slopes)
Q10_vals <- summary(slopes)
Q10_vals$Q10 <- exp(10 * Q10_vals$Temperature.trend)

#also calculate Q10 for subsets: for the basalt dataset
dataCO2FB_C <- data_original[which(data_original$Treatment=="C1"|data_original$Treatment=="C3"|data_original$Treatment=="FB3"|data_original$Treatment=="FB1"),]
dataCO2FB_C$logCO2 <- log(dataCO2FB_C$µmolCO2m2s)
model_Q10 <- lmer(logCO2 ~ Temperature * rock + (1 | PotNo),
                  data = dataCO2FB_C)
summary(model_Q10)
Anova(model_Q10)

slopes <- emtrends(model_Q10, ~ rock, var = "Temperature")
# Pairwise tests: Are Q10 values different?
pairs(slopes)
Q10_vals <- summary(slopes)
Q10_vals$Q10 <- exp(10 * Q10_vals$Temperature.trend)




# # If Temperature:rock interaction is significant → Q10 differs among rocks
# # Extract slopes per rock (b parameter)
# slopes <- emtrends(model_Q10, ~ rock, var = "Temperature")
# # Pairwise tests: Are Q10 values different?
# pairs(slopes)
# Q10_vals <- summary(slopes)
# Q10_vals$Q10 <- exp(10 * Q10_vals$Temperature.trend)

#3% SOC soils only
dataCO2D_B_size <- data_original[which(data_original$SOC==3),] 
dataCO2D_B_size$logCO2 <- log(dataCO2D_B_size$µmolCO2m2s)
model_Q10 <- lmer(logCO2 ~ Temperature * rock + (1 | PotNo),
                  data = dataCO2D_B_size)
summary(model_Q10)
Anova(model_Q10)

slopes <- emtrends(model_Q10, ~ rock, var = "Temperature")
# Pairwise tests: Are Q10 values different?
pairs(slopes)
Q10_vals <- summary(slopes)
Q10_vals$Q10 <- exp(10 * Q10_vals$Temperature.trend)




```



```{r}
data_original$SWC2 = data_original$SWC^2
data_original_complete <- data_original[complete.cases (data_original$SWC),]
data_original_complete <- data_original_complete[complete.cases (data_original_complete$temperature),]
data_original_complete <- data_original_complete[complete.cases (data_original_complete$µmolCO2m2s),]


results <- data.frame(
  PotNo = numeric(0),
  Treatment = character(0),
  b = numeric(0),
  SE_b = numeric(0)
)

for (p in unique(data_original_complete$PotNo)) {
  
  # Subset data for one pot
  subset_data <- subset(data_original_complete, PotNo == p)

  # Log-linear respiration model
  modelCO2 <- lm(log(µmolCO2m2s) ~ Temperature + SWC + SWC2, data = subset_data)

  # Extract slope and SE
  b     <- coef(modelCO2)["Temperature"]
  SE_b  <- summary(modelCO2)$coefficients["Temperature", "Std. Error"]

  treatment <- unique(subset_data$Treatment)

  # Store
  results <- rbind(results, data.frame(
    PotNo = p,
    Treatment = treatment,
    b = b,
    SE_b = SE_b
  ))
}

results$Q10     <- exp(10 * results$b)
results$SE_Q10  <- results$Q10 * 10 * results$SE_b

library(Rmisc)

Q10_treat_summary <- summarySE(
  data = results,
  measurevar = "Q10",
  groupvars = "Treatment",
  na.rm = TRUE
)

Q10_treat_summary$CI_low  <- Q10_treat_summary$Q10 - Q10_treat_summary$ci
Q10_treat_summary$CI_high <- Q10_treat_summary$Q10 + Q10_treat_summary$ci




# results <- data.frame(
#   PotNo = numeric(0),
#   b = numeric(0),
#   Q10 = numeric(0),
#   Q10_low = numeric(0),
#   Q10_high = numeric(0)
# )
# 
# for (p in unique(data_original_complete$PotNo)) {
#   
#   # Subset data for this pot
#   subset_data <- subset(data_original_complete, PotNo == p)
#   # Fit model on the subset
#   
#   # 
#   modelCO2 <- lm(log(µmolCO2m2s) ~ Temperature+ SWC + SWC2 , data = subset_data )
# 
#   summary (modelCO2)
#   # Extract the temperature coefficient
#  # Extract b (temperature slope)
#   b <- coef(modelCO2)["Temperature"]
#   
#   # Q10
#   Q10 <- exp(10 * b)
#   # CI for b
#   ci_b <- confint(modelCO2)["Temperature", ]
#   b_low <- ci_b[1]
#   b_high <- ci_b[2]
#   
#   treatment <- unique (subset_data$Treatment)
#   
#   # Compute Q10
#   Q10 <- exp(10 * b)
#   
#   # Confidence interval for b
#   ci_b <- confint(modelCO2, parm = "Temperature", method = "Wald")
#   b_low <- ci_b[1]
#   b_high <- ci_b[2]
#   
#   # Convert CI of b to CI of Q10
#   Q10_low <- exp(10 * b_low)
#   Q10_high <- exp(10 * b_high)
#   
#   # Add to results dataframe
#   results <- rbind(results, data.frame(
#     PotNo = p,
#     b = b,
#     Q10 = Q10,
#     Q10_low = Q10_low,
#     Q10_high = Q10_high,
#     Treatment = treatment
#   ))
# }
# 
# 
# library(Rmisc)
# results_summary <- summarySE(
#   data = results,
#   measurevar = "Q10",
#   groupvars = "Treatment"
# )
# 
# # for (p in unique(data_original_complete$PotNo)) {
# #   subset_data <- subset(data_original_complete, PotNo == p)
# #   modelCO2 <- lmer (data = data_original, µmolCO2m2s ~ Temperature + SWC + SWC2  + (1 | PotNo) )
# # summary (modelCO2)
# # coef (modelCO2)
# # b <- fixef(modelCO2)["Temperature"]
# # Q10 = exp (10*b)
# #  
# # 
# # 
# # }
# 
library(ggplot2)
# 
Q10_treat_summary1 <- Q10_treat_summary[which(Q10_treat_summary$Treatment=="C1"|Q10_treat_summary$Treatment=="FB1" ),]
Q10_treat_summary3 <- Q10_treat_summary[-which(Q10_treat_summary$Treatment=="C1"|Q10_treat_summary$Treatment=="FB1" ),]
# 
Q10plot_1 <- ggplot(Q10_treat_summary1 , aes(x = Treatment, y = Q10)) +
  geom_point(size = 3) +
  geom_errorbar(aes(ymin = CI_low, ymax = CI_high), width = 0.2) +
  labs(
    x = "Treatment",
    y = "Q10 (± CI)",
  ) +
  theme_minimal(base_size = 14)


Q10plot_3 <- ggplot(Q10_treat_summary3 , aes(x = Treatment, y = Q10)) +
  geom_point(size = 3) +
  geom_errorbar(aes(ymin = CI_low, ymax = CI_high), width = 0.2) +
  labs(
    x = "Treatment",
    y = "Q10 (± CI)",
  ) +
  theme_minimal(base_size = 14)

library(cowplot)

Q10_subplot <- plot_grid(
  Q10plot_1, 
  Q10plot_3,
  labels = c("A", "B"),
  ncol = 2,
  align = "hv"
)


ggsave(
  filename = "Q10_subplot.jpeg",
  plot = Q10_subplot,
  width = 12,
  height = 6,
  dpi = 600
)




```

