##-------------------------------------------------## ## Introduction to Structural Equation Modelling ## ## John Fox ## ## FIOCRUZ/PROCC Brasil ## ## November 2008 ## ##-------------------------------------------------## library(sem) # Duncan, Haller, and Portes's nonrecursive peer-influences model # FIML estimates from the sem function R.DHP <- read.moments(diag=FALSE, names=c('ROccAsp', 'REdAsp', 'FOccAsp', 'FEdAsp', 'RParAsp', 'RIQ', 'RSES', 'FSES', 'FIQ', 'FParAsp')) .6247 .3269 .3669 .4216 .3275 .6404 .2137 .2742 .1124 .0839 .4105 .4043 .2903 .2598 .1839 .3240 .4047 .3054 .2786 .0489 .2220 .2930 .2407 .4105 .3607 .0186 .1861 .2707 .2995 .2863 .5191 .5007 .0782 .3355 .2302 .2950 .0760 .0702 .2784 .1988 .1147 .1021 .0931 -.0438 .2087 # Note: lower triangle of correlation matrix # path parameter start-value model.DHP.1 <- specify.model() RIQ -> ROccAsp, gamma51, NA RSES -> ROccAsp, gamma52, NA FSES -> FOccAsp, gamma63, NA FIQ -> FOccAsp, gamma64, NA FOccAsp -> ROccAsp, beta56, NA ROccAsp -> FOccAsp, beta65, NA ROccAsp <-> ROccAsp, sigma77, NA FOccAsp <-> FOccAsp, sigma88, NA ROccAsp <-> FOccAsp, sigma78, NA model.DHP.1 sem.DHP.1 <- sem(model.DHP.1, R.DHP, 329, fixed.x=c('RIQ', 'RSES', 'FSES', 'FIQ')) summary(sem.DHP.1) # Blau and Duncan's recursive basic-stratification model R.BD <- read.moments(diag=FALSE, names=c("FatherEd", "FatherOcc", "Education", "FirstJob", "1962Job")) .516 .453 .438 .332 .417 .538 .322 .405 .596 .541 R.BD model.BD <- specify.model() FatherEd -> Education, gamma31, NA FatherOcc -> Education, gamma32, NA FatherOcc -> FirstJob, gamma41, NA Education -> FirstJob, beta43, NA FatherOcc -> 1962Job, gamma52, NA Education -> 1962Job, beta53, NA FirstJob -> 1962Job, beta54, NA Education <-> Education, sigma66, NA FirstJob <-> FirstJob, sigma77, NA 1962Job <-> 1962Job, sigma88, NA model.BD sem.BD <- sem(model.BD, R.BD, 20700, fixed.x=c("FatherEd", "FatherOcc")) summary(sem.BD) # note: the OLS estimates # Duncan, Haller, and Portes's block-recursive peer-influences model model.DHP.2 <- specify.model() RIQ -> ROccAsp, gamma51, NA RSES -> ROccAsp, gamma52, NA FSES -> ROccAsp, gamma53, NA RSES -> FOccAsp, gamma62, NA FSES -> FOccAsp, gamma63, NA FIQ -> FOccAsp, gamma64, NA FOccAsp -> ROccAsp, beta56, NA ROccAsp -> FOccAsp, beta65, NA RIQ -> REdAsp, gamma71, NA RSES -> REdAsp, gamma72, NA FSES -> REdAsp, gamma73, NA RSES -> FEdAsp, gamma82, NA FSES -> FEdAsp, gamma83, NA FIQ -> FEdAsp, gamma84, NA ROccAsp -> REdAsp, beta75, NA FOccAsp -> FEdAsp, beta86, NA FEdAsp -> REdAsp, beta78, NA REdAsp -> FEdAsp, beta87, NA ROccAsp <-> ROccAsp, sigma99, NA FOccAsp <-> FOccAsp, sigma1010, NA ROccAsp <-> FOccAsp, sigma910, NA REdAsp <-> REdAsp, sigma1111, NA FEdAsp <-> FEdAsp, sigma1212, NA REdAsp <-> FEdAsp, sigma1112, NA model.DHP.2 sem.DHP.2 <- sem(model.DHP.2, R.DHP, 329, fixed.x=c('RIQ', 'RSES', 'FSES', 'FIQ')) summary(sem.DHP.2) # Duncan, Haller and Portes peer-influences model with multiple indicators # of latent endogenous variables R.DHP <- read.moments(diag=FALSE, names=c('ROccAsp', 'REdAsp', 'FOccAsp', 'FEdAsp', 'RParAsp', 'RIQ', 'RSES', 'FSES', 'FIQ', 'FParAsp')) .6247 .3269 .3669 .4216 .3275 .6404 .2137 .2742 .1124 .0839 .4105 .4043 .2903 .2598 .1839 .3240 .4047 .3054 .2786 .0489 .2220 .2930 .2407 .4105 .3607 .0186 .1861 .2707 .2995 .2863 .5191 .5007 .0782 .3355 .2302 .2950 .0760 .0702 .2784 .1988 .1147 .1021 .0931 -.0438 .2087 # Note: lower triangle of correlation matrix (repeated) model.DHP.3 <- 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 # fixed parameter RGenAsp -> REdAsp, lam21, NA FGenAsp -> FOccAsp, NA, 1 # fixed parameter 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 model.DHP.3 sem.DHP.3 <- sem(model.DHP.3, R.DHP, 329, fixed.x=c('RParAsp', 'RIQ', 'RSES', 'FSES', 'FIQ', 'FParAsp')) summary(sem.DHP.3) # the model with equality constraints on corresponding parameters model.DHP.4 <- specify.model() 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, beta, NA, RGenAsp -> FGenAsp, beta, NA, RGenAsp -> ROccAsp, NA, 1, # fixed parameter RGenAsp -> REdAsp, lam, NA, FGenAsp -> FOccAsp, NA, 1, # fixed parameter FGenAsp -> FEdAsp, lam, NA, RGenAsp <-> RGenAsp, ps11, NA, FGenAsp <-> FGenAsp, ps11, NA, RGenAsp <-> FGenAsp, ps12, NA, ROccAsp <-> ROccAsp, theta1, NA, REdAsp <-> REdAsp, theta2, NA, FOccAsp <-> FOccAsp, theta1, NA, FEdAsp <-> FEdAsp, theta2, NA model.DHP.4 sem.DHP.4 <- sem(model.DHP.4, R.DHP, 329, fixed.x=c('RParAsp', 'RIQ', 'RSES', 'FSES', 'FIQ', 'FParAsp')) summary(sem.DHP.4) # likelihood-ratio test for equality constraints anova(sem.DHP.3, sem.DHP.4) # Confirmatory Factor Analyis of Holzinger's data on 9 psychological tests R.Holzinger <- read.moments(diag=FALSE, names=c("Word.meaning", "Sentence.completion", "Odd.words", "Mixed.arithmetic", "Remainders", "Missing.numbers", "Gloves", "Boots", "Hatchets")) .75 .78 .72 .44 .52 .47 .45 .53 .48 .82 .51 .58 .54 .82 .74 .21 .23 .28 .33 .37 .35 .30 .32 .37 .33 .36 .38 .45 .31 .30 .37 .31 .36 .38 .52 .67 R.Holzinger model.Holzinger.1 <- specify.model() # model assuming uncorrelated factors Verbal -> Word.meaning, lam11, NA Verbal -> Sentence.completion, lam21, NA Verbal -> Odd.words, lam31, NA Arithmetic -> Mixed.arithmetic, lam42, NA Arithmetic -> Remainders, lam52, NA Arithmetic -> Missing.numbers, lam62, NA Spatial -> Gloves, lam73, NA Spatial -> Boots, lam83, NA Spatial -> Hatchets, lam93, NA Word.meaning <-> Word.meaning, th1, NA Sentence.completion <-> Sentence.completion,th2, NA Odd.words <-> Odd.words, th3, NA Mixed.arithmetic <-> Mixed.arithmetic, th4, NA Remainders <-> Remainders, th5, NA Missing.numbers <-> Missing.numbers, th6, NA Gloves <-> Gloves, th7, NA Boots <-> Boots, th8, NA Hatchets <-> Hatchets, th9, NA Verbal <-> Verbal, NA, 1 Arithmetic <-> Arithmetic, NA, 1 Spatial <-> Spatial, NA, 1 model.Holzinger.1 sem.Holzinger.1 <- sem(model.Holzinger.1, R.Holzinger, 696) summary(sem.Holzinger.1) model.Holzinger.2 <- specify.model() # model allowing correlated factors Verbal -> Word.meaning, lam11, NA Verbal -> Sentence.completion, lam21, NA Verbal -> Odd.words, lam31, NA Arithmetic -> Mixed.arithmetic, lam42, NA Arithmetic -> Remainders, lam52, NA Arithmetic -> Missing.numbers, lam62, NA Spatial -> Gloves, lam73, NA Spatial -> Boots, lam83, NA Spatial -> Hatchets, lam93, NA Word.meaning <-> Word.meaning, th1, NA Sentence.completion <-> Sentence.completion,th2, NA Odd.words <-> Odd.words, th3, NA Mixed.arithmetic <-> Mixed.arithmetic, th4, NA Remainders <-> Remainders, th5, NA Missing.numbers <-> Missing.numbers, th6, NA Gloves <-> Gloves, th7, NA Boots <-> Boots, th8, NA Hatchets <-> Hatchets, th9, NA Verbal <-> Verbal, NA, 1 Arithmetic <-> Arithmetic, NA, 1 Spatial <-> Spatial, NA, 1 Verbal <-> Arithmetic, psy12, NA Verbal <-> Spatial, psy13, NA Arithmetic <-> Spatial, psy23, NA sem.Holzinger.2 <- sem(model.Holzinger.2, R.Holzinger, 696) summary(sem.Holzinger.2) # test of factor correlations anova(sem.Holzinger.1, sem.Holzinger.2)