##### SEM.5.2 - ANIMAL EXAMPLE OF MEDIATION - DROUGHT EFFECTS ON RODENT INTERACTIONS

library(lavaan)
library(AICcmodavg)


# load dataset		
sem.dat=read.csv("SEM_5_2_Animal Example of Mediation_data.csv")
dim(sem.dat)
summary(sem.dat)

	# Each row provides rodent count data at a given plot in a given year. Data are pooled for this illustration.
	
	# DIndex is Drought Index 1: 7 months of weather data were aggregated into this index.
	
	# 'Rare' and 'Giant': counts were ln(...+1) transformed (i.e. using in R log())
	
	# Predictors (DIndex and Giant) were lag values from previous years. Note that when year=2007, rodent data are missing as surveys started in 2006.


######## Models ##########
## Net Effect Model (added by jbg)
  # net effect
    net <- 'Rare ~ DIndex'
    # Fit model
    net.fit <- sem(net, data=sem.dat)
    # Output results
    summary(net.fit, standardized=T,rsq=T)

## Test of mediation 

	# Complete mediation
		complete <- 'Giant ~ DIndex
				         Rare ~ Giant'
		# Fit model
		complete.fit <- sem(complete, data=sem.dat)
		# Output results
		summary(complete.fit,standardized=T,rsq=T)

	# Partial mediation
		partial <- 'Giant ~ DIndex
				        Rare ~  Giant + DIndex'
		partial.fit <- sem(partial, data=sem.dat)
		summary(partial.fit,standardized=T,rsq=T)

	# No mediation
		none <- 'Giant ~ 0*DIndex					
				     Rare ~  0*Giant + DIndex'
		none.fit <- sem(none, data=sem.dat)
		summary(none.fit, rsq=T, standardized=T)	

# Test of Mediation:
	anova(complete.fit,partial.fit,none.fit)	# partial mediation model is best, as shown in Fig 3.b.
	
# AICc comparison table
	aictab(c(complete.fit,partial.fit,none.fit), c("Complete", "Partial", "No Mediation"))
	

## Indirect and Total Effects
	# Partial mediation
	partial <- 'Giant ~ a*DIndex
	            Rare ~  b*Giant + c*DIndex

	            direct   := c          
              indirect := a*b
              total    := c + (a*b)
            '
	
	
	partial.fit <- sem(partial, data=sem.dat)
	summary(partial.fit,standardized=T,rsq=T)