#2022-12-28 Modified script to also output SEM evaluation as .tsv files. /André Jernung, Swedish National Data Service
#Load packages
library(nlme)
library(piecewiseSEM)

#Set working directory to location of script
setwd(".")

#Define dataset
df1 <- read.delim("DataWeatherVegMoose.tsv")

#Create SEM overview table
SEM_overview <- data.frame(matrix(ncol=10,nrow=0, dimnames=list(NULL, c("SEM.object.name", "included.LMEs", "Fisher.C", "df", "P.Value", "AIC", "AICc", "BIC", "K", "n"))))

#Run piecewise SEM. Proportion of days over 20 degrees C and fireweed
modelEPA_20 <- psem(
  lme(NDF_Epa ~ PropDaysOver20 + TotalPrecipGS, random = ~ 1 | Site, correlation=corCAR1(form=~Year|Site),na.action = na.omit,data = df1),
  lme(N_Epa ~ PropDaysOver20 + TotalPrecipGS, random = ~ 1 | Site, correlation=corCAR1(form=~Year|Site),na.action = na.omit,data = df1),
  lme(MeanWt ~ PropDaysOver20 + TotalPrecipGS + N_Epa + NDF_Epa, random = ~ 1 | Site, correlation=corCAR1(form=~Year|Site),na.action = na.omit,data = df1),
  NDF_Epa %~~% N_Epa
)
summary(modelEPA_20, standardize = "scale")

#Output coefficients and R2 tables
modelEPA_20_out <- summary(modelEPA_20, standardize = "scale")
colnames(modelEPA_20_out$coefficients)[9] <- "Signif.codes"
write.table(modelEPA_20_out["coefficients"], file='modelEPA_20_coefficients.tsv', quote=FALSE, sep='\t', col.names=colnames(modelEPA_20_out$coefficients), row.names=FALSE)
write.table(modelEPA_20_out["R2"], file='modelEPA_20_R2.tsv', quote=FALSE, sep='\t', col.names=colnames(modelEPA_20_out$R2), row.names=FALSE)

#Add LMEs and statistics to SEM overview table
SEM_LMEs <- gsub("\n ", " |", modelEPA_20_out$call)
SEM_overview[nrow(SEM_overview) + 1,] <- c(modelEPA_20_out$name, SEM_LMEs, modelEPA_20_out$Cstat$Fisher.C, modelEPA_20_out$Cstat$df, modelEPA_20_out$Cstat$P.Value, modelEPA_20_out$IC$AIC, modelEPA_20_out$IC$AICc, modelEPA_20_out$IC$BIC, modelEPA_20_out$IC$K, modelEPA_20_out$IC$n)

#Run piecewise SEM. Mean average daily temperature and fireweed
modelEPA_mean <- psem(
  lme(NDF_Epa ~ MeanAvgTempGS + TotalPrecipGS, random = ~ 1 | Site, correlation=corCAR1(form=~Year|Site),na.action = na.omit,data = df1),
  lme(N_Epa ~ MeanAvgTempGS + TotalPrecipGS, random = ~ 1 | Site, correlation=corCAR1(form=~Year|Site),na.action = na.omit,data = df1),
  
  lme(MeanWt ~ MeanAvgTempGS + TotalPrecipGS + N_Epa + NDF_Epa, random = ~ 1 | Site, correlation=corCAR1(form=~Year|Site),na.action = na.omit,data = df1),
  NDF_Epa %~~% N_Epa
)
summary(modelEPA_mean, standardize = "scale")

#Output coefficients and R2 tables
modelEPA_mean_out <- summary(modelEPA_mean, standardize = "scale")
colnames(modelEPA_mean_out$coefficients)[9] <- "Signif.codes" #This column does not have a heading in the data.frame
write.table(modelEPA_mean_out["coefficients"], file='modelEPA_mean_coefficients.tsv', quote=FALSE, sep='\t', col.names=colnames(modelEPA_mean_out$coefficients), row.names=FALSE)
write.table(modelEPA_mean_out["R2"], file='modelEPA_mean_R2.tsv', quote=FALSE, sep='\t', col.names=colnames(modelEPA_mean_out$R2), row.names=FALSE)

#Add LMEs and statistics to SEM overview table
SEM_LMEs <- gsub("\n ", " |", modelEPA_mean_out$call)
SEM_overview[nrow(SEM_overview) + 1,] <- c(modelEPA_mean_out$name, SEM_LMEs, modelEPA_mean_out$Cstat$Fisher.C, modelEPA_mean_out$Cstat$df, modelEPA_mean_out$Cstat$P.Value, modelEPA_mean_out$IC$AIC, modelEPA_mean_out$IC$AICc, modelEPA_mean_out$IC$BIC, modelEPA_mean_out$IC$K, modelEPA_mean_out$IC$n)


#Run piecewise SEM. Proportion of days over 20 degrees C and downy birch
modelBPU_20 <- psem(
  lme(NDF_Bpu ~ PropDaysOver20 + TotalPrecipGS, random = ~ 1 | Site, correlation=corCAR1(form=~Year|Site),na.action = na.omit,data = df1),
  lme(N_Bpu ~ PropDaysOver20 + TotalPrecipGS, random = ~ 1 | Site, correlation=corCAR1(form=~Year|Site),na.action = na.omit,data = df1),
  
  lme(MeanWt ~ PropDaysOver20 + TotalPrecipGS + N_Bpu + NDF_Bpu, random = ~ 1 | Site, correlation=corCAR1(form=~Year|Site),na.action = na.omit,data = df1),
  NDF_Bpu %~~% N_Bpu
)
summary(modelBPU_20, standardize = "scale")

#Output coefficients and R2 tables
modelBPU_20_out <- summary(modelBPU_20, standardize = "scale")
colnames(modelBPU_20_out$coefficients)[9] <- "Signif.codes" #This column does not have a heading in the data.frame
write.table(modelBPU_20_out["coefficients"], file='modelBPU_20_coefficients.tsv', quote=FALSE, sep='\t', col.names=colnames(modelBPU_20_out$coefficients), row.names=FALSE)
write.table(modelBPU_20_out["R2"], file='modelBPU_20_R2.tsv', quote=FALSE, sep='\t', col.names=colnames(modelBPU_20_out$R2), row.names=FALSE)

#Add LMEs and statistics to SEM overview table
SEM_LMEs <- gsub("\n ", " |", modelBPU_20_out$call)
SEM_overview[nrow(SEM_overview) + 1,] <- c(modelBPU_20_out$name, SEM_LMEs, modelBPU_20_out$Cstat$Fisher.C, modelBPU_20_out$Cstat$df, modelBPU_20_out$Cstat$P.Value, modelBPU_20_out$IC$AIC, modelBPU_20_out$IC$AICc, modelBPU_20_out$IC$BIC, modelBPU_20_out$IC$K, modelBPU_20_out$IC$n)

#Run piecewise SEM. Mean average daily temperature and downy birch
modelBPU_mean <- psem(
  lme(NDF_Bpu ~ MeanAvgTempGS + TotalPrecipGS, random = ~ 1 | Site, correlation=corCAR1(form=~Year|Site),na.action = na.omit,data = df1),
  lme(N_Bpu ~ MeanAvgTempGS + TotalPrecipGS, random = ~ 1 | Site, correlation=corCAR1(form=~Year|Site),na.action = na.omit,data = df1),
  lme(MeanWt ~ MeanAvgTempGS + TotalPrecipGS + N_Bpu + NDF_Bpu, random = ~ 1 | Site, correlation=corCAR1(form=~Year|Site),na.action = na.omit,data = df1),
  NDF_Bpu %~~% N_Bpu
)
summary(modelBPU_mean, standardize = "scale")

#Output coefficients and R2 tables
modelBPU_mean_out <- summary(modelBPU_mean, standardize = "scale")
colnames(modelBPU_mean_out$coefficients)[9] <- "Signif.codes" #This column does not have a heading in the data.frame
write.table(modelBPU_mean_out["coefficients"], file='modelBPU_mean_coefficients.tsv', quote=FALSE, sep='\t', col.names=colnames(modelBPU_mean_out$coefficients), row.names=FALSE)
write.table(modelBPU_mean_out["R2"], file='modelBPU_mean_R2.tsv', quote=FALSE, sep='\t', col.names=colnames(modelBPU_mean_out$R2), row.names=FALSE)

#Add LMEs and statistics to SEM overview table
SEM_LMEs <- gsub("\n ", " |", modelBPU_mean_out$call)
SEM_overview[nrow(SEM_overview) + 1,] <- c(modelBPU_mean_out$name, SEM_LMEs, modelBPU_mean_out$Cstat$Fisher.C, modelBPU_mean_out$Cstat$df, modelBPU_mean_out$Cstat$P.Value, modelBPU_mean_out$IC$AIC, modelBPU_mean_out$IC$AICc, modelBPU_mean_out$IC$BIC, modelBPU_mean_out$IC$K, modelBPU_mean_out$IC$n)

#Output SEM overview table
write.table(SEM_overview, file='SEM_overview.tsv', quote=FALSE, sep='\t', col.names=colnames(SEM_overview), row.names=FALSE)

cat("\n\n")

#Output R environment for logging purposes
sessionInfo()