library(nlme)
library(lme4)
library(piecewiseSEM)
library(QuantPsyc)

karstdata <- read.csv("data.csv",header = T )

# Standardize the data
karst_data_matrix <- scale(karstdata)

# Convert the matrix to a data frame
Rt <- as.data.frame(karst_data_matrix)

# Keep the column names from the original data
colnames(Rt) <- colnames(karstdata)

### 1. Climate
# composite the climate variables
model1 <- lm(CUE ~  MAP + MAT, Rt)
coefs(model1, standardize = 'scale')
beta_MAP <-  summary(model1)$coefficients[2, 1]
beta_MAT <- summary(model1)$coefficients[3, 1]

climate <- beta_MAP * Rt$MAP + beta_MAT * Rt$MAT
Rt$climate <- climate

### 2. Microbial properties
# composite the microbial variables
model2 <- lm(CUE ~  Shannon_16S + Abundance_16S + Kr_16S +
              Shannon_ITS + Abundance_ITS + Kr_ITS, Rt)
coefs(model2, standardize = 'scale')

beta_Shannon_16S <-  summary(model2)$coefficients[2, 1]
beta_Abundance_16S <- summary(model2)$coefficients[3, 1]
beta_Kr_16S <- summary(model2)$coefficients[4, 1]
beta_Shannon_ITS <- summary(model2)$coefficients[5, 1]
beta_Abundance_ITS <- summary(model2)$coefficients[6, 1]
beta_Kr_ITS <- summary(model2)$coefficients[7, 1]

microbe <- beta_Shannon_16S * Rt$Shannon_16S + beta_Abundance_16S * Rt$Abundance_16S +  beta_Kr_16S * Rt$Kr_16S +
          beta_Shannon_ITS * Rt$Shannon_ITS + beta_Abundance_ITS * Rt$Abundance_ITS +  beta_Kr_ITS * Rt$Kr_ITS
Rt$microbe <- microbe


### 3. soil minerals
#composite the mineral variables
model3 <- lm(CUE ~  Caexe + CaCO3 + Fed, Rt)
coefs(model3, standardize = 'scale')

beta_Caexe <-  summary(model3)$coefficients[2, 1]
beta_CaCO3 <-  summary(model3)$coefficients[3, 1]
beta_Fed <-  summary(model3)$coefficients[4, 1]

minerals <- beta_Caexe * Rt$Caexe + beta_CaCO3 * Rt$CaCO3 +  beta_Fed * Rt$Fed
Rt$minerals <- minerals


## piecewise SEM
CUE.list2 <- list(
  lm(CUE ~ microbe + minerals + Macro , na.action = na.omit, data = Rt),
  lm(microbe ~ minerals + Macro, na.action = na.omit, data = Rt),
  lm(Macro ~ AGB + minerals + climate, na.action = na.omit, data = Rt),
  lm(minerals ~ AGB + climate, na.action = na.omit, data = Rt),
  lm(AGB ~ climate, na.action = na.omit, data = Rt)
)
CUE.psem2 <- as.psem(CUE.list2)
(new.summary2 <- summary(CUE.psem2, .progressBar = F))
plot(CUE.psem2)