
#Necessary packages
library(lme4) 
library(car)
library(mgcv)
library(ggplot2)
library(moments)


#only with dataset completed for temperature, soil water content (SWC) and  Flux (with greater than 0 values)
GHGData <- GHGData[which(GHGData$Flux > 0),] 


#remove outliers based on the histogram to correct the data skewness without a data transformation
hist(GHGData$Flux)
hist(GHGData$temperature)
hist(GHGData$SWC)

GHGData <- GHGData[which(GHGData$Flux > n),]
GHGData <- GHGData[which(GHGData$Flux < n),]
GHGData <- GHGData[which(GHGData$SWC > n),]

#check the skewness of  flux, and SWC data
skewness(GHGData$Flux, na.rm=TRUE) 
skewness(GHGData$SWC, na.rm=TRUE) 


##Model with lmer - normalize the effects of temperature and SWC on the changes in flux

#models with log10 data transformation
model1 = lmer(log10(Flux) ~ temperature + SWC + (1|Mesocosm_no),data=GHGData)
model2 = lmer(log10(Flux) ~ temperature + log10(SWC + 1) + (1|Mesocosm_no),data=GHGData)
model3 = lmer(log10(Flux) ~ temperature + log10(SWC^2 + SWC +1) + (1|Mesocosm_no),data=GHGData)
model4 = lmer(log10(Flux) ~ temperature + SWC + SWC^2 + (1|Mesocosm_no),data=GHGData)

#models without data transformation
model5 = lmer(Flux ~ temperature + SWC + (1|Mesocosm_no),data=GHGData)
model6 = lmer(Flux ~ temperature + log10(SWC + 1) + (1|Mesocosm_no),data=GHGData)
model7 = lmer(Flux ~ temperature + log10(SWC^2 + SWC +1) + (1|Mesocosm_no),data=GHGData)
model8 = lmer(Flux ~ temperature + SWC + SWC^2 + (1|Mesocosm_no),data=GHGData)

#check for low AIC values to find out good models

AIC(model1,model2,model3,model4,model5,model6,model7,model8)

summary(MODEL)

#check the normality of the model visually (with qqnorm and hist functions) and mathematically (with Shapiro test)
qqnorm(resid(MODEL))
hist(resid(MODEL))
shapiro.test(residuals(MODEL))

##Model with lmer - Find out the effect of treatments (Basalt and B. subtilis) and time (Day) on the temperature and SWC normalized  flux data


GHGData$residualsSCE <- residuals (MODEL)

model = lmer(residualsSCE ~  Day*Basalt_bin*Bacillus_bin + (1|Mesocosm_no),data=GHGData)


summary(model)
vif(model)
shapiro.test(residuals(model))

#statistical analysis with ANOVA
Anova(model, type="II")




