A first look at structured equation models using the Lavaan package

SEM example.R

# A FIRST LOOK AT LAVAAN
## By Grace Charles, presented at Davis R Users' Group on May 15, 2015
## adapted from Jim Grace's SEM workshop and Lavaan tutorials

# Set your working directory
setwd("~/Desktop/DAVIS/sem workshop")


###Load Libraries
library(lavaan)
library(qgraph)

#Built in dataset

data(PoliticalDemocracy)
PD<-PoliticalDemocracy
head(PD)

## here is an example model from one of Lavaan's build-in datasets
## the measurement model equations are "latent" and represented by =~
## regressions are indicated by ~
## residual correlations (in this case because they represent different years of the same measurement) are represented by ~~

model <- '
  # measurement model
    ind60 =~ x1 + x2 + x3
    dem60 =~ y1 + y2 + y3 + y4
    dem65 =~ y5 + y6 + y7 + y8
  # regressions
    dem60 ~ ind60
    dem65 ~ ind60 + dem60
  # residual correlations
    y1 ~~ y5
    y2 ~~ y4 
    y2~~ y6
    y3 ~~ y7
    y4 ~~ y8
    y6 ~~ y8
'

#fit your SEM
fit <- sem(model, data = PD)

#summarize results
summary(fit, standardized = TRUE, rsq = T)

##plot results using semPaths function in qgraph
semPaths(fit, "std", edge.label.cex = 0.5, curvePivot = TRUE, layout = "tree")

##check to see if you missed anything. High mi values suggest that there is a path that you missed.

modindices(fit)

## looks good

##can also look at variance tables  
vartable(fit)


## sometimes you get warnings about the scale of your variables
#Warning message:
# In getDataFull(data = data, group = group, group.label = group.label, :
#                 lavaan WARNING: some observed variances are (at least) a factor 100 times larger than others; please rescale

# in this case, all you have to do to make this error go away is rescale variables

#model comparison

#you can compare alternative pathway models using AIC, BIC, etc:

#create second alternative model
names(PD)
  
model2 <- '
  # measurement model
    ind60 =~ x1 + x2 + x3
    dem60 =~ y1 + y2 + y3 + y4
    dem65 =~ y5 + y6 + y7 + y8
  # regressions
    dem60 ~ ind60
    dem65 ~ dem60
#took out ind60 from regression
  # residual correlations
    y1 ~~ y5
    y2 ~~ y4 + y6
    y3 ~~ y7
    y4 ~~ y8
    y6 ~~ y8
'
fit2 <- sem(model2, data = PD)
summary(fit2)

AIC(fit, fit2)

## what about nonlinear data?

# Set your working directory
("~/Desktop/sem workshop")
#commands in bold

# Load data and name file ?k.dat?
k.dat<-read.csv("./Keeley_rawdata_select4.csv")

# Examine contents of keeley data file
names(k.dat)
head(k.dat)


# Write lavaan code for this single equation model
mod <- '
    rich ~ cover
    cover ~ firesev
'

k.dat$cov2<-k.dat$cover^2

mod2<- '
      rich ~ cover + cov2
      cover ~ firesev
      cover ~~ cov2
      cov2 ~~ firesev
'



# Fit the model (i.e. est. parameters)
mod1.fit <- sem(mod, data=k.dat)
mod2.fit<- sem(mod2, data=k.dat,fixed.x=FALSE)

#need to rescale data.
vartable(mod1.fit)
k.dat$rich<-k.dat$rich/100

# Output a summary of the computed results - summary of mod2 suggests that both cover and cover squared can impact 
summary(mod1.fit, rsq=T)  # rsq=T means output the r-sqr
summary(mod2.fit, rsq=T)


semPaths(mod1.fit, "std", edge.label.cex = 0.5, curvePivot = TRUE, layout = "tree")
semPaths(mod2.fit, "std", edge.label.cex = 0.5, curvePivot = TRUE, layout = "tree")

add

library(semPlot)

It happens to me too when I try to load the semPlot package. Is the problem solved?

Thanks, but where is the "Keeley_rawdata_select4.csv" file. Could we use this file too please?