####### SEM.6-SEM VERSUS MULTIPLE REGRESSION v2
####### code to match "SEM_6_SEM_versus_Multiple_Regression_v2.pptx"
####### revised 18.08.17

### Load data, recode data
k.dat <- read.csv("SEM_6_SEM_versus_Multiple_Regression_data.csv")
names(k.dat)

### Recode vars to roughly same scale
k6.dat <- subset(k.dat, select=c(rich, firesev, abiotic, hetero, age, distance))
summary(k6.dat); cov(k6.dat)
k6.dat$rich     <- k6.dat$rich/10
k6.dat$firesev  <- k6.dat$firesev
k6.dat$abiotic  <- k6.dat$abiotic/10
k6.dat$hetero   <- k6.dat$hetero*10
k6.dat$age      <- k6.dat$age/10
k6.dat$distance <- k6.dat$distance/10


##### Multiple regression models using lm
## Initial model including all predictors
mr.lm1 <- lm(rich ~ firesev + abiotic + hetero + age + distance, data=k6.dat); summary(mr.lm1)

## Simplified by removing age
mr.lm2 <- lm(rich ~ firesev + abiotic + hetero + distance, data=k6.dat); summary(mr.lm2)

## Compare models with and without age as predictor
# classic anova 
anova(mr.lm1, mr.lm2) # results indicate no significant loss of information removing age

# AICc comparison
library(AICcmodavg)
aictab(list(mr.lm1, mr.lm2), c("lm_model1", "lm_model2")) # results indicate model 2 slightly better


########### SEM FOR FIRE RECOVERY STUDY #########
library(lavaan)
## First run most comprehensive model "mod.2" 
# and check for missing paths

# specify "mod.2" 
mod.2 <- 'rich ~ abiotic + hetero + distance 
                 + firesev + age
          abiotic ~ distance
          hetero ~ distance
          age ~ distance
          firesev ~ age'

# Estimate model "mod.2"
mod.2.fit <- sem(mod.2, data=k6.dat)
summary(mod.2.fit)


## Compare to model omitting age -> rich, i.e., "mod.1" 
mod.1 <- 'rich ~ abiotic + hetero + distance 
                 + firesev
          abiotic ~ distance
          hetero ~ distance
          age ~ distance
          firesev ~ age'

# Estimate model "mod.1"
mod.1.fit <- sem(mod.1, data=k6.dat)


### Model comparison based on AIC
# Compare models using anova function
anova(mod.1.fit, mod.2.fit)

# Compare models using AICc table
aictab(list(mod.1.fit, mod.2.fit), 
       c("mod.1", "mod.2"))

# Detailed results for model 1
summary(mod.1.fit)


### Request R-squares for model1
summary(mod.1.fit, rsq=T)

### Request std coefs for model1
summary(mod.1.fit, standardized=T)
standardizedsolution(mod.1.fit)