# Code for analysis of soil properties and vegetation cover change in precipitation experiment using R program

#Load packages required for mixed models and estimating p-vales from mixed models

library(lme4)
library(languageR)

## Test normality of soil variables
 
hist(organics)  
hist(clay)
hist(sand)
hist(silt)

shapiro.test(organics) 
shapiro.test(clay)
shapiro.test(sand)
shapiro.test(silt)

### Blocked Anova for comparisons of soil properties between cheatgrass and wheatgrass dominated soils
## Blocked anovas also performed for soil carbon, soil nitrogen, and CN ratios using the same code
organic.anova <- aov(organics ~ soilid + block)
clay.anova <- aov(clay ~ soilid + block)
sand.anova <- aov(sand ~ soilid + block)
silt.anova <- aov(silt ~ soilid + block)

###Linear mixed models (for NO3, NH4, and total soil nitrogen) in soils over time

nitrogen_model = lmer(NO3 ~ soilid + (1|date) + (1|block), data=data1)
summary(nitrogen_model)


####### Linear mixed models for vegetation cover variables ####

# Variables that should be factors are read as factors

ftreat <- factor(veg_june$treatment)
fyear <- factor(veg_june$year)


### make a variable for random effect of plot within block

fplot <- factor(veg_june$plot, labels="p")
fblock <- factor(veg_june$block, labels="b")
plotblock <- rep(NA, nrow(veg_june))

for (i in 1:nrow(veg_june)){
  plotblock[i] <- paste(fplot[i], fblock[i], sep="")
}

####Check response variables, tranform cover values that violate assumptions of normality ####
## Introduced forbs, Native forbs, and warm-season grasses were all square-root transformed using the code below ##

hist(native_forb)
sqnative <- sqrt(native_forb)
hist(sqnative)
shapiro.test(sqnative)

# Re-organize years so 2014 is compared to all other years

years <- factor(fyear, levels=c("2014","2011","2012","2013"))

## subset data to look only at treatment differences in 2014

only14 <- veg_june[which(year=="2014"),] 

### Below are linear mixed-effect models for exotic grasses in June. Introduced forbs, warm-season grasses, native forbs,
### and cool-season grasses in both June and AUgust were analysed with models using the same structure

####### Exotic grasses - lmer for treatment and time ######## 

xgrass <- lmer(ex_grass ~ ftreat + years + (1|fblock) + (1|plotblock), data=veg_june)
summary(xgrass)
pvals.fnc(xgrass)


##### Exotic grasses - lmer for comparison of treatments in 2014 only ####

xxgrass <- lmer(ex_grass ~ treatment + (1|block), data=only14)
summary(xxgrass)
pvals.fnc(xxgrass)


#### model diagnostics, run for all models to check assumptions of normality of variances #####

fit <-fitted(xgrass)
r <- resid(xgrass)


par(mfrow=c(2,2))

#histogram of residuals 

hist(r, xlab = "Residuals", main = "Histogram of residuals",freq=FALSE)
rr <- seq(min(r),max(r),length.out=100)
lines(rr,dnorm(rr,mean=0,sd=sd(r)),col="red")
box()

#residuals vs. fitted

plot(fit,r, ylab = "Residuals", xlab = "Fitted values",
     main = "Residuals vs fitted")
abline(h=0,col="red")


#quantile-quantile plot
qqnorm(r)
qqline(r)




### paired t tests for nitrogen addition experiment
#data first checked for normality and transformed if not normal
# checked to see whether variances were equal
> var.test(totnat_a, totnat_n)

#paired t test- when variances were equal 
t.test(totnat_a, totnat_n, paired = TRUE, var.equal=TRUE)

#paired  t test- when variances were unequal 
t.test(totnat_a, totnat_n, paired = TRUE)