R14

		Structural Equation Models

library(sem)	# Requires installation

JS.cor <- as.matrix(read.table("H:/S5600/Data sets/ASCII/Ch10/Job_Sat.txt"))
JS.cor


# CFA - Measures reliable for Job sat, self esteem?

model.JS.cfa <- specify.model()		# Set up your model
JS -> Job_sat1,		lam11, NA	# Common Factor loadings
JS -> Job_sat2,		lam21, NA
SE -> Self_est1,	lam32, NA
SE -> Self_est2,	lam42, NA
Job_sat1 <-> Job_sat1,	th1, NA		# Specific Factor Variances
Job_sat2 <-> Job_sat2,	th2, NA
Self_est1 <-> Self_est1,th3, NA
Self_est2 <-> Self_est2,th4, NA
JS <-> SE,		ph12, NA	# Common Factor Covariance(s)
JS <-> JS,		NA, 1		# Common Factor Variances (fixed to 1)
SE <-> SE,		NA, 1
					# Blank line to end input

model.JS.cfa
n <- 106

cfa.JS <- sem(model.JS.cfa, JS.cor, n)

cfa.JS
summary(cfa.JS)

# SEM - Job sat as a function of self esteem

model.JS.sem <- specify.model()		# Set up your model
JS -> Job_sat1,		NA, 1		# Factor loadings
JS -> Job_sat2,		lamy2, NA	# Note Job_sat1 fixed to 1
SE -> Self_est1,	lamx1, NA
SE -> Self_est2,	lamx2, NA
SE -> JS,		gam, NA	
Job_sat1 <-> Job_sat1,	eps1, NA	# Specific Factor Variances
Job_sat2 <-> Job_sat2,	eps2, NA
Self_est1 <-> Self_est1,del1, NA
Self_est2 <-> Self_est2,del2, NA
JS <-> JS,		psi, NA		# Common Factor Variances
SE <-> SE,		NA, 1		# JS Variance not fixed (endogenous)
					# Blank line to end input

model.JS.sem

sem.JS <- sem(model.JS.sem, JS.cor, n)

sem.JS
summary(sem.JS)


# standardized solution

sem.JS$coeff
v <- sem.JS$coeff[4]^2+sem.JS$coeff[9]	# variance of JS - scale to 1
v

gam.s <- sem.JS$coeff[4]/sqrt(v)
gam.s
psi.s <- sem.JS$coeff[9]/v
psi.s

gam.s^2 + psi.s		# should equal 1 now

lamy1.s <- sqrt(v)	# standardized lambdas
lamy1.s
lamy2.s <- sem.JS$coeff[1]*sqrt(v)
lamy2.s

std.coef(sem.JS)


# New data - Duncan, Haller, Portes data on student and friend 
# aspirations

DHP.cor <- as.matrix(read.table("F:/S5600/Data sets/DHP.txt"))
DHP.cor

model.DHP.sem1 <- specify.model()
         RParAsp  -> RGenAsp, gam11,  NA
         RIQ      -> RGenAsp, gam12,  NA
         RSES     -> RGenAsp, gam13,  NA
         FSES     -> RGenAsp, gam14,  NA
         RSES     -> FGenAsp, gam23,  NA
         FSES     -> FGenAsp, gam24,  NA
         FIQ      -> FGenAsp, gam25,  NA
         FParAsp  -> FGenAsp, gam26,  NA
         FGenAsp  -> RGenAsp, beta12, NA
         RGenAsp  -> FGenAsp, beta21, NA
         RGenAsp  -> ROccAsp,  NA,     1	# One loading fixed to 1
         RGenAsp  -> REdAsp,  lam21,  NA	# For each factor
         FGenAsp  -> FOccAsp,  NA,     1
         FGenAsp  -> FEdAsp,  lam42,  NA
         RGenAsp <-> RGenAsp, ps11,   NA
         FGenAsp <-> FGenAsp, ps22,   NA
         RGenAsp <-> FGenAsp, ps12,   NA
         ROccAsp <-> ROccAsp, theta1, NA
         REdAsp  <-> REdAsp,  theta2, NA
         FOccAsp <-> FOccAsp, theta3, NA
         FEdAsp  <-> FEdAsp,  theta4, NA


sem1.DHP <- sem(model.DHP.sem1, DHP.cor, 329,  
         fixed.x=c('RParAsp', 'RIQ', 'RSES', 'FSES', 'FIQ', 'FParAsp'))

# fixed.x -> purely input (exogenous) variables; saves us from having
# to specify their variances in the model

summary(sem1.DHP)

# Drop gam14, gam23 loadings (nonsignificant)

model.DHP.sem2 <- specify.model()
         RParAsp  -> RGenAsp, gam11,  NA
         RIQ      -> RGenAsp, gam12,  NA
         RSES     -> RGenAsp, gam13,  NA
         FSES     -> FGenAsp, gam24,  NA
         FIQ      -> FGenAsp, gam25,  NA
         FParAsp  -> FGenAsp, gam26,  NA
         FGenAsp  -> RGenAsp, beta12, NA
         RGenAsp  -> FGenAsp, beta21, NA
         RGenAsp  -> ROccAsp,  NA,     1
         RGenAsp  -> REdAsp,  lam21,  NA
         FGenAsp  -> FOccAsp,  NA,     1
         FGenAsp  -> FEdAsp,  lam42,  NA
         RGenAsp <-> RGenAsp, ps11,   NA
         FGenAsp <-> FGenAsp, ps22,   NA
         RGenAsp <-> FGenAsp, ps12,   NA
         ROccAsp <-> ROccAsp, theta1, NA
         REdAsp  <-> REdAsp,  theta2, NA
         FOccAsp <-> FOccAsp, theta3, NA
         FEdAsp  <-> FEdAsp,  theta4, NA

sem2.DHP <- sem(model.DHP.sem2, DHP.cor, 329,  
         fixed.x=c('RParAsp', 'RIQ', 'RSES', 'FSES', 'FIQ', 'FParAsp'))

summary(sem2.DHP)

# Test differences (nested model)

pchisq(29.997 - 26.697, 17 - 15, lower.tail=F)

# Reduce model by forcing equality of equivalent loadings

model.DHP.sem3 <- specify.model()
         RParAsp  -> RGenAsp, gam1,  NA
         RIQ      -> RGenAsp, gam2,  NA
         RSES     -> RGenAsp, gam3,  NA
         FSES     -> FGenAsp, gam3,  NA
         FIQ      -> FGenAsp, gam2,  NA
         FParAsp  -> FGenAsp, gam1,  NA
         FGenAsp  -> RGenAsp, beta, NA
         RGenAsp  -> FGenAsp, beta, NA
         RGenAsp  -> ROccAsp,  NA,     1
         RGenAsp  -> REdAsp,  lam,  NA
         FGenAsp  -> FOccAsp,  NA,     1
         FGenAsp  -> FEdAsp,  lam,  NA
         RGenAsp <-> RGenAsp, psi,   NA
         FGenAsp <-> FGenAsp, psi,   NA
         RGenAsp <-> FGenAsp, ps12,   NA
         ROccAsp <-> ROccAsp, theta1, NA
         REdAsp  <-> REdAsp,  theta2, NA
         FOccAsp <-> FOccAsp, theta1, NA
         FEdAsp  <-> FEdAsp,  theta2, NA

sem3.DHP <- sem(model.DHP.sem3, DHP.cor, 329,  
         fixed.x=c('RParAsp', 'RIQ', 'RSES', 'FSES', 'FIQ', 'FParAsp'))

summary(sem3.DHP)

pchisq(36.529 - 29.997, 25 - 17, lower.tail=F)

mod.indices(sem3.DHP)

# Add correlation between respondent and friend's 
# occupational aspirations

model.DHP.sem4 <- specify.model()
         RParAsp  -> RGenAsp, gam1,  NA
         RIQ      -> RGenAsp, gam2,  NA
         RSES     -> RGenAsp, gam3,  NA
         FSES     -> FGenAsp, gam3,  NA
         FIQ      -> FGenAsp, gam2,  NA
         FParAsp  -> FGenAsp, gam1,  NA
         FGenAsp  -> RGenAsp, beta, NA
         RGenAsp  -> FGenAsp, beta, NA
         RGenAsp  -> ROccAsp,  NA,     1
         RGenAsp  -> REdAsp,  lam,  NA
         FGenAsp  -> FOccAsp,  NA,     1
         FGenAsp  -> FEdAsp,  lam,  NA
         RGenAsp <-> RGenAsp, psi,   NA
         FGenAsp <-> FGenAsp, psi,   NA
         RGenAsp <-> FGenAsp, ps12,   NA
         ROccAsp <-> ROccAsp, theta1, NA
         REdAsp  <-> REdAsp,  theta2, NA
         FOccAsp <-> FOccAsp, theta1, NA
         FEdAsp  <-> FEdAsp,  theta2, NA
	 FOccAsp <-> ROccAsp, theta24, NA
	 
sem4.DHP <- sem(model.DHP.sem4, DHP.cor, 329,  
         fixed.x=c('RParAsp', 'RIQ', 'RSES', 'FSES', 'FIQ', 'FParAsp'))

summary(sem4.DHP)

pchisq(36.529 - 26.193, 25 - 24, lower.tail=F)

std.coef(sem4.DHP)