install.packages(c("MASS","sem"),contriburl = "http://cran.stat.ucla.edu/bin/windows/contrib/2.6/")
library(MASS)
library(sem)


#Path Analysis 1: Simple mediation model
cordata <- read.moments(diag=FALSE, names=c('FamilyBackground', 'Ability', 'Motivation', 'Coursework', 'Achievement'))
     .417
     .190  .205
     .372  .498  .375
     .417  .737  .255   .615

#For the following, arrows will indicate paths, followed by an arbitrary path name.  NA tells the program to estimate starting values for the coefficients.
model1 <- specify.model()
    FamilyBackground	-> Achievement, gam1, NA
    FamilyBackground	-> Ability, gam2, NA
    Ability		-> Achievement, gam3, NA
    FamilyBackground <-> FamilyBackground, psi1, NA
    Ability <-> Ability, psi2, NA
    Achievement <-> Achievement, psi3, NA

path.model1 <- sem(model1, cordata, 1000)  #will give a warning just to note that we aren't using all variables available

summary(path.model1)

#Compare coefficients with regression output.  The first part creates simulated standardized data with the exact same correlations.
cormat=matrix(c(1,.417,.190,.372,.417,.417,1,.205,.498,.737,.190,.205,1,.375,.255,
.372,.498,.375,1,.615,.417,.737,.255,.615,1),5, byrow=T)
cormat  

newdata= data.frame(mvrnorm(n=1000,mu=c(0,0,0,0,0),Sigma=cormat, empirical=T))
names(newdata)=c("FamilyBackground", "Ability", "Motivation", "Coursework", "Achievement")

(pathgam2=lm(Ability~FamilyBackground, newdata))
(pathgam23=lm(Achievement~Ability+FamilyBackground, newdata))
summary(path.model1)

#Path Model 2: A more complex model. Add motivation as a second mediator of the Family Background - Achievement relationship.

model2 <- specify.model()
    FamilyBackground	-> Achievement, gam11, NA
    FamilyBackground	-> Ability, gam12, NA
    FamilyBackground	-> Motivation, gam13, NA
    Ability		-> Achievement, gam21, NA
    Motivation    -> Achievement, gam31, NA
    FamilyBackground <-> FamilyBackground, psi1, NA
    Ability <-> Ability, psi2, NA
    Motivation <-> Motivation, psi3, NA
    Achievement <-> Achievement,  psi4, NA


path.model2 <- sem(model2, cordata, 1000)  #To bypass 'estimating' FamilyBackground, could put fixed.x="FamilyBackground" in this function and take its variance (psi) out of the above model
summary(path.model2)  #As this model is overidentified, we now have several fit indices available to us


#Path Model 3: To path model 2, add a third mediator Coursework, however, it makes sense that Ability and Motivation would have an effect on that as well.
model3 <- specify.model()
    FamilyBackground	-> Achievement, gam11, NA
    FamilyBackground	-> Ability, gam12, NA
    FamilyBackground	-> Motivation, gam13, NA
    FamilyBackground	-> Coursework, gam14, NA
	Ability		-> Achievement, gam21, NA
	Ability		-> Coursework,  gam22, NA
	Motivation        -> Coursework,  gam31, NA
	Motivation        -> Achievement,  gam32, NA
    	Coursework        -> Achievement, gam41, NA 
    FamilyBackground <-> FamilyBackground, psi1, NA
    Ability <-> Ability, psi2, NA
    Motivation <-> Motivation, psi3, NA
    Coursework <-> Coursework, psi4, NA
    Achievement <-> Achievement,  psi5, NA

path.model3 <- sem(model3, cordata, 1000)  

summary(path.model3)      

#Path Model 4: A final model for comparison.  Paths added based on modification indices.  Leaves only an indirect path for Motivation - Achievement and suggests Ability causes Motivation
model4 <- specify.model()
    FamilyBackground	-> Achievement, gam11, NA
    FamilyBackground	-> Ability, gam12, NA
    FamilyBackground	-> Motivation, gam13, NA
    FamilyBackground	-> Coursework, gam14, NA
	Ability		-> Achievement, gam21, NA
	Ability		-> Coursework,  gam22, NA
	Ability		-> Motivation,  gam23, NA		
	Motivation        -> Coursework,  gam31, NA
    	Coursework        -> Achievement, gam41, NA 
    FamilyBackground <-> FamilyBackground, psi1, NA
    Ability <-> Ability, psi2, NA
    Motivation <-> Motivation, psi3, NA
    Coursework <-> Coursework, psi4, NA
    Achievement <-> Achievement,  psi5, NA

path.model4 <- sem(model4, cordata, 1000)  

summary(path.model4)     #Best fit of all.