################## SEM.2.2-Local Estimation with piecewiseSEM Rscript #########
# Load data - using same data file as in SEM_2_1_Intro to Lavaan
dat <- read.csv("SEM.2.1_data.csv")
attach(dat)

## Recode vars so results should match with SEM_2_1_Intro to Lavaan
x1 <- x1/100
y1 <- y1/100
y2 <- y2
y3 <- y3/100
t.dat <- data.frame(x1, y1, y2, y3)


##### Demonstrating Principles by Hand

# Specify model as series of regressions
y1.mod <- lm(y1 ~ x1, data=t.dat)
y2.mod <- lm(y2 ~ x1, data=t.dat)
y3.mod <- lm(y3 ~ y1 + y2, data=t.dat)

# Capture residuals
x1.res <- x1
y1.res <- resid(y1.mod)
y2.res <- resid(y2.mod)
y3.res <- resid(y3.mod)

# Look at conditional independences
dev.new(height=4, width=8); par(mfrow=c(1,2))
plot(y3.res ~ x1.res, pch=16); cor.test(y3.res, x1.res); abline(lm(y3.res ~ x1.res), lwd=2)
plot(y2.res ~ y1.res, pch=16); cor.test(y2.res, y1.res)

# Respecify model
y1.mod <- lm(y1 ~ x1, data=t.dat)
y2.mod <- lm(y2 ~ x1, data=t.dat)
y3.mod2 <- lm(y3 ~ y1 + y2 + x1, data=t.dat)


##### Using Lefcheck's piecewiseSEM procedures
library(piecewiseSEM)

# Model 1 - SE model is a list of the local models
pw.mod1 = psem(
  lm(y1 ~ x1, t.dat),
  lm(y2 ~ x1, t.dat),
  lm(y3 ~ y1 + y2, t.dat))

# request full summary
summary(pw.mod1, t.dat)

# request just the dSep tests
dSep(pw.mod1)


# Revised Model based on dSep test
pw.mod2 = psem(
  lm(y1 ~ x1, t.dat),
  lm(y2 ~ x1, t.dat),
  lm(y3 ~ y1 + y2 + x1, t.dat))

# request just the dSep tests
dSep(pw.mod2)
fisherC(pw.mod2)

# request full summary
summary(pw.mod2)

# request coefficients
coefs(pw.mod2, standardize="none")

# compare models
AIC(pw.mod1, pw.mod2, aicc=T)