################################################
# Path models
################################################
rm(list = ls())
model_data <- read_excel("model_data.xlsx")
model_data$CoverShrubLayer_trans <- as.vector(scale(asin(sqrt(model_data$CoverShrubLayer))))
model_data$TSF_trans <- 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$herb_q0_scaled <- as.numeric(scale(model_data$herb_q0))
model_data$pH_scaled <- as.numeric(scale(model_data$pH))

## Species richness ----
model_data_sem <- model_data[, c("CoverShrubLayer_trans",
                                 "TSF_trans",
                                 "canopy_num",
                                 "fence_num",
                                 "pH_scaled",
                                 "plot",
                                 "herb_q0_scaled")]

any(is.na(model_data_sem))

### Shrubs
sem.shrub<- lmer(CoverShrubLayer_trans ~
                   canopy_num +
                   fence_num +
                   pH_scaled +
                   (1|plot),
                 data = model_data_sem)
### TSF
sem.TSF<- lmer(TSF_trans ~ 
                 canopy_num +
                 CoverShrubLayer_trans + 
                 (1|plot),
               data = model_data_sem)

# Species richness
sem.SN<- lmer(herb_q0_scaled ~ 
                TSF_trans +
                CoverShrubLayer_trans +
                fence_num +
                pH_scaled +
                (1|plot), 
              data = model_data_sem)

## Combining models
model = psem(
  sem.shrub,
  sem.TSF,
  sem.SN)
summary(model)

# SES MPD ----
names(model_data)
model_data_sem <- model_data[, c("CoverShrubLayer_trans",
                                 "TSF_trans",
                                 "canopy_num",
                                 "fence_num",
                                 "pH_scaled",
                                 "plot",
                                 "ses_scaled")]

any(is.na(model_data_sem))
sum(is.na(model_data_sem)) # 1 plots have to be removed
model_data_sem <- na.omit(model_data_sem)

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

### TSF
sem.TSF<- lmer(TSF_trans ~ 
                 canopy_num +
                 CoverShrubLayer_trans + 
                 (1|plot),
               data = model_data_sem)

### SES MPD
sem.SESMPD <- lmer(ses_scaled ~ 
                     TSF_trans +
                     fence_num +
                     CoverShrubLayer_trans +
                     pH_scaled +
                     (1|plot), 
                   data = model_data_sem)

## Combining models
model = psem(
  sem.shrub,
  sem.TSF,
  sem.SESMPD)
summary(model)