##------------------------------------------------------------------------## ## Script for web appendix on structural equation models ## ## An R and S-PLUS Companion to Applied Regression ## ## John Fox ## ## Sage Publications, 2002 ## ##------------------------------------------------------------------------## # preliminaries options(width=65) options(digits=5) # two-stage least squares library(sem) data(Klein) Klein attach(Klein) p.lag <- c(NA, p[-length(p)]) x.lag <- c(NA, x[-length(x)]) a <- year - 1931 cbind(year, a, p, p.lag, x, x.lag) eqn.1 <- tsls(c ~ p + p.lag + I(wp + wg), instruments= ~ g + t + wg + a + p.lag + k.lag + x.lag) summary(eqn.1) eqn.2 <- tsls(i ~ p + p.lag + k.lag, instruments= ~ g + t + wg + a + p.lag + k.lag + x.lag) summary(eqn.2) eqn.3 <- tsls(wp ~ x + x.lag + a, instruments= ~ g + t + wg + a + p.lag + k.lag + x.lag) summary(eqn.3) # general structural-equation model R.dhp <- matrix(c( # lower triangle of correlation matrix 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, .6247, 1, 0, 0, 0, 0, 0, 0, 0, 0, .3269, .3669, 1, 0, 0, 0, 0, 0, 0, 0, .4216, .3275, .6404, 1, 0, 0, 0, 0, 0, 0, .2137, .2742, .1124, .0839, 1, 0, 0, 0, 0, 0, .4105, .4043, .2903, .2598, .1839, 1, 0, 0, 0, 0, .3240, .4047, .3054, .2786, .0489, .2220, 1, 0, 0, 0, .2930, .2407, .4105, .3607, .0186, .1861, .2707, 1, 0, 0, .2995, .2863, .5191, .5007, .0782, .3355, .2302, .2950, 1, 0, .0760, .0702, .2784, .1988, .1147, .1021, .0931, -.0438, .2087, 1 ), ncol=10, byrow=T) rownames(R.dhp) <- colnames(R.dhp) <- c('ROccAsp', 'REdAsp', 'FOccAsp', 'FEdAsp', 'RParAsp', 'RIQ', 'RSES', 'FSES', 'FIQ', 'FParAsp') R.dhp ram.dhp <- matrix(c( # arrow parameter start-value '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', 'bet12', NA, 'RGenAsp -> FGenAsp', 'bet21', NA, 'RGenAsp -> ROccAsp', NA, 1, 'RGenAsp -> REdAsp', 'lamy21', NA, 'FGenAsp -> FOccAsp', NA, 1, 'FGenAsp -> FEdAsp', 'lamy42', NA, 'RGenAsp <-> RGenAsp', 'psi11', NA, 'FGenAsp <-> FGenAsp', 'psi22', NA, 'RGenAsp <-> FGenAsp', 'psi12', NA, 'ROccAsp <-> ROccAsp', 'theps1', NA, 'REdAsp <-> REdAsp', 'theps2', NA, 'FOccAsp <-> FOccAsp', 'theps3', NA, 'FEdAsp <-> FEdAsp', 'theps4', NA), ncol=3, byrow=T) ram.dhp sem.dhp <- sem(ram.dhp, R.dhp, N=329, fixed.x=c('RParAsp', 'RIQ', 'RSES', 'FSES', 'FIQ', 'FParAsp')) sem.dhp summary(sem.dhp) # a model with equality constraints ram.dhp.1 <- matrix(c( 'RParAsp -> RGenAsp', 'gam1', NA, 'RIQ -> RGenAsp', 'gam2', NA, 'RSES -> RGenAsp', 'gam3', NA, 'FSES -> RGenAsp', 'gam4', NA, 'RSES -> FGenAsp', 'gam4', NA, 'FSES -> FGenAsp', 'gam3', NA, 'FIQ -> FGenAsp', 'gam2', NA, 'FParAsp -> FGenAsp', 'gam1', NA, 'FGenAsp -> RGenAsp', 'bet', NA, 'RGenAsp -> FGenAsp', 'bet', NA, 'RGenAsp -> ROccAsp', NA, 1, 'RGenAsp -> REdAsp', 'lamy', NA, 'FGenAsp -> FOccAsp', NA, 1, 'FGenAsp -> FEdAsp', 'lamy', NA, 'RGenAsp <-> RGenAsp', 'psi', NA, 'FGenAsp <-> FGenAsp', 'psi', NA, 'RGenAsp <-> FGenAsp', 'psi12', NA, 'ROccAsp <-> ROccAsp', 'theps1', NA, 'REdAsp <-> REdAsp', 'theps2', NA, 'FOccAsp <-> FOccAsp', 'theps1', NA, 'FEdAsp <-> FEdAsp', 'theps2', NA), ncol=3, byrow=T) sem.dhp.1 <- sem(ram.dhp.1, R.dhp, N=329, fixed.x=c('RParAsp', 'RIQ', 'RSES', 'FSES', 'FIQ', 'FParAsp')) summary(sem.dhp.1) mod.indices(sem.dhp.1) ram.dhp.2 <- matrix(c( 'RParAsp -> RGenAsp', 'gam1', NA, 'RIQ -> RGenAsp', 'gam2', NA, 'RSES -> RGenAsp', 'gam3', NA, 'FSES -> RGenAsp', 'gam4', NA, 'RSES -> FGenAsp', 'gam4', NA, 'FSES -> FGenAsp', 'gam3', NA, 'FIQ -> FGenAsp', 'gam2', NA, 'FParAsp -> FGenAsp', 'gam1', NA, 'FGenAsp -> RGenAsp', 'bet', NA, 'RGenAsp -> FGenAsp', 'bet', NA, 'RGenAsp -> ROccAsp', NA, 1, 'RGenAsp -> REdAsp', 'lamy', NA, 'FGenAsp -> FOccAsp', NA, 1, 'FGenAsp -> FEdAsp', 'lamy', NA, 'RGenAsp <-> RGenAsp', 'psi', NA, 'FGenAsp <-> FGenAsp', 'psi', NA, 'RGenAsp <-> FGenAsp', 'psi12', NA, 'ROccAsp <-> ROccAsp', 'theps1', NA, 'REdAsp <-> REdAsp', 'theps2', NA, 'FOccAsp <-> FOccAsp', 'theps1', NA, 'FEdAsp <-> FEdAsp', 'theps2', NA, 'FOccAsp <-> ROccAsp', 'theps24', NA), ncol=3, byrow=T) sem.dhp.2 <- sem(ram.dhp.2, R.dhp, N=329, fixed.x=c('RParAsp', 'RIQ', 'RSES', 'FSES', 'FIQ', 'FParAsp')) summary(sem.dhp.2)