####### SEM.5-THE TEST OF MEDIATION

### Set working directory, load data, recode data
setwd("E:/ppt_files/_education/Code-Data-Materials") #commands in bold
k.dat <- read.csv("SEM_5_The Test of Mediation_data_Keeley_rawdata_select4.csv")
names(k.dat)

### Recode vars to roughly same scale
k3.dat <- subset(k.dat, select=c(cover, firesev, age))
mean(k3.dat[,1]); mean(k3.dat[,2]); mean(k3.dat[,3]); cov(k3.dat)
k3.dat$cover <- k3.dat$cover*10
k3.dat$firesev <- k3.dat$firesev
k3.dat$age <- k3.dat$age/10

### Load Libraries
library(lavaan)
library(AICcmodavg)


######### TEST OF MEDIATION ##########
### Net (total) effect of age on cover
mod.1 <- 'cover ~ age'                 # Specify model
mod.1.fit <- sem(mod.1, data=k3.dat)   # Fit the model
summary(mod.1.fit, stand=T, rsq=T)     # Extract results

### Three degrees of mediation
## Complete mediation
comp.mod  <- 'cover ~ firesev
              firesev ~ age'
comp.mod.fit <- sem(comp.mod, data=k3.dat)

## Partial mediation
partial.mod <- 'cover ~ firesev + age
                firesev ~ age'
partial.mod.fit <- sem(partial.mod, data=k3.dat)

## No mediation
nomed.mod <- 'firesev ~ 0*age
              cover ~ age + 0*firesev'
nomed.mod.fit <- sem(nomed.mod, data=k3.dat)

## Compare models with anova
anova(comp.mod.fit, partial.mod.fit, nomed.mod.fit)

## Compare three models with AICc
aictab.lavaan(list(comp.mod.fit, partial.mod.fit, nomed.mod.fit), c("Complete", "Partial", "None"))

## Compare two models with AICc
aictab.lavaan(list(comp.mod.fit, partial.mod.fit), c("Complete", "Partial"))

## Complete mediation alternative specification - does not seem to be useful
comp.mod2  <- 'cover ~ firesev + 0*age
               firesev ~ age'
comp.mod2.fit <- sem(comp.mod2, data=k3.dat)

## Compare two models with AICc  - same result
aictab.lavaan(list(comp.mod2.fit, partial.mod.fit), c("Complete2", "Partial"))

# a small digression: asking for the intercepts
comp.mod.fit <- sem(comp.mod, meanstructure=T, data=k3.dat)
summary(comp.mod.fit)

### Compute indirect and total effects for partial mediation model
## Partial mediation
mod.4 <- 'cover ~ b*firesev + c*age
          firesev ~ a*age

          # define quantities 
          direct   := c 
          indirect := a*b
          total    := c + (a*b)
         '

# Fit the model
mod.4.fit <- sem(mod.4, data=k.dat)

# Extract results
summary(mod.4.fit, stand=T, rsq=T)