### Code for conducting SEM in R
### Uses package 'sem'

library(sem)
help(sem)

### import data
datum=read.csv(file.choose())
head(datum)

### first step is to specify the path model, which is done with 'specify.model'
help(specifyModel)

###Build the path diagram
###Arguments are path, variable name, starting value (NA lets R choose starting value)
modelPop=specifyModel()
	CWD<>CWD, CWDVar, NA ### variation in course woody debris
	SmMDens<>SmMDens, MVar, NA ### error in small mammal density
	WeasDens<>WeasDens, WVar, NA ### error in weasel density
	CWD>WeasDens, WeasVeg, NA ### direct path between CWD and weasel density
	CWD>SmMDens, PreyVeg, NA ### indirect path through CWD and weasel density 
	SmMDens>WeasDens, WeasPrey, NA ### through effects of CWD on small mammal density
# Note that model is fully specified:
### # parameters = # observations


### Run the SEM

results=sem(modelPop,data=datum)


### summary provides goodness-of-fit tests
### parameter estimates for variances/errors are standard deviation estimates of residuals
### parameter estimates for paths are betas of coefficients in regression equations
### compare estiamtes to truth (in excel data-sheet)
### NOTE - THESE ARE NOT STANDARDIZED PARAMETER COEFFICIENTS 
summary(results)


### function std.coef provides standardized path coefficiets
help(stdCoef)
stdCoef(results)
### Path coefficients for variances/errors are % of data not explained by path model (1-r^2)
### Path coefficients for paths are correlation coefficients (r) after accounting for other effects

### compare to correlation matrix
cor(datum)

#Note that PreyVeg = r for CWD>SmMDens
#Note that PreyVeg = Sqrt of R-square for Endogenous variables for SmMDens
#Note that WeasVeg+PreyVeg*WesPrey = r for CWD>WeasDens
### Why is R^2 for WeasDens greater than square of r for CWD>WeasDens?
### Because WeasDens is also a function of SmMamDens, which is a function of more than just CWD



### run a simpler model without effect of CWD on small mammals
modelPop2=specifyModel()
	CWD<>CWD, CWDVar, NA
	SmMDens<>SmMDens, MVar, NA
	WeasDens<>WeasDens, WVar, NA
	CWD>SmMDens, PreyVeg, NA
	SmMDens>WeasDens, WeasPrey, NA
#

results2=sem(modelPop2,data=datum)
summary(results2)

###Note additional goodness of fit measures