Structural Equation Modeling in R

谢益辉 2007-04-23

I feel unhappy, so I just decide to write something not-so-interesting.

Everybody is using Structural Equation Models (SEM) in China as if it were shameful that you didn’t construct SEMs in your paper. Again, few people know SEM theories – not to say the computation issues. They don’t care, because they know this is a so-called advanced tool. This morning I’ve talked about several ironic things about SEM in my speech.

SEM Paths

The theory can be regarded as easy, though in fact it’s not. The main idea is to develop a model with a theoretical covariance matrix ($\Sigma$) which is similar to the sample covariance matrix S. Therefore, how to minimize their discrepancy becomes the most critical part in parameter estimation. Common discrepancy functions can be derived through maximum likelihood (ML), unweighted least squares (ULS) and generalized least squares (GLS), etc. Take ML as an example, the discrepancy function is

$$FML = tr(S\Sigma-1) + \log|\Sigma| - log|S| - (p + q)$$

This function is quite complicated; usually it has more than 20 parameters. All along I have been doubting the reliability of its optimization (the dimension is really too high!). However, nobody ever mentioned this problem. Maybe they’re too happy to notice it.

OK, since my topic is the implementation in R, let’s check the package “sem” first. The function sem() is used to fit general structural equation models.

sem {sem} R Documentation 

**General Structural Equation Models **

**Description** 

`sem` fits general structural equation models (with both observed and unobserved variables) by the method of maximum likelihood, assuming multinormal errors. Observed variables are also called _indicators_ or _manifest variables_; unobserved variables are also called _factors_ or _latent variables_. Normally, the generic function (`sem`) would be used. 


**Usage**
    
    sem(ram, ...)
    
    ## S3 method for class 'mod':
    sem(ram, S, N, obs.variables=rownames(S), fixed.x=NULL, debug=FALSE, ...)
        
    ## Default S3 method:
    sem(ram, S, N, param.names = paste("Param", 1:t, sep = ""), 
        var.names = paste("V", 1:m, sep = ""), fixed.x = NULL, raw=FALSE, 
        debug = FALSE, analytic.gradient = TRUE, warn = FALSE, maxiter = 500, 
        par.size=c('ones', 'startvalues'), refit=TRUE, start.tol=1E-6, ...) 
        
    startvalues(S, ram, debug = FALSE, tol=1E-6)
    
    ## S3 method for class 'sem':
    print(x, ...)
    
    ## S3 method for class 'sem':
    summary(object, digits=5, conf.level=0.9, ...)
    
    ## S3 method for class 'sem':
    deviance(object, ...)
    
    ## S3 method for class 'sem':
    df.residual(object, ...)

Currently only ML and_ two-stage least squares (TSLS) _are available in this package. The specification is easy: just tell R your sample covariance matrix and settings for coefficients and covariances. The examples are quite clear:

# Note: These examples can't be run via example() because the default file
#  argument of specify.model() requires that the model specification be entered
#  at the command prompt. The examples can be copied and run in the R console,
#  however. See ?specify.model for further information

    ## Not run: 

# ------------- Duncan, Haller and Portes peer-influences model ----------------------
# A nonrecursive SEM with unobserved endogenous variables and fixed exogenous variables

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=TRUE)
            
rownames(R.DHP) <- colnames(R.DHP) <- 
    c('ROccAsp', 'REdAsp', 'FOccAsp', 'FEdAsp', 'RParAsp', 'RIQ',
        'RSES', 'FSES', 'FIQ', 'FParAsp')

# Fit the model using a symbolic ram specification

model.dhp <- 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
    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

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

summary(sem.dhp.1)

##     Model Chisquare =  26.697   Df =  15 Pr(>Chisq) = 0.031302
##     Chisquare (null model) =  872   Df =  45
##     Goodness-of-fit index =  0.98439
##     Adjusted goodness-of-fit index =  0.94275
##     RMSEA index =  0.048759   90 
##     Bentler-Bonnett NFI =  0.96938
##     Tucker-Lewis NNFI =  0.95757
##     Bentler CFI =  0.98586
##     BIC =  -60.244 
##    
##     Normalized Residuals
##       Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##    -0.7990 -0.1180  0.0000 -0.0120  0.0397  1.5700 
##    
##     Parameter Estimates
##           Estimate  Std Error z value Pr(>|z|)                       
##    gam11   0.161224 0.038487   4.1890 2.8019e-05 RGenAsp <--- RParAsp
##    gam12   0.249653 0.044580   5.6001 2.1428e-08 RGenAsp <--- RIQ    
##    gam13   0.218404 0.043476   5.0235 5.0730e-07 RGenAsp <--- RSES   
##    gam14   0.071843 0.050335   1.4273 1.5350e-01 RGenAsp <--- FSES   
##    gam23   0.061894 0.051738   1.1963 2.3158e-01 FGenAsp <--- RSES   
##    gam24   0.228868 0.044495   5.1437 2.6938e-07 FGenAsp <--- FSES   
##    gam25   0.349039 0.044551   7.8346 4.6629e-15 FGenAsp <--- FIQ    
##    gam26   0.159535 0.040129   3.9755 7.0224e-05 FGenAsp <--- FParAsp
##    beta12  0.184226 0.096207   1.9149 5.5506e-02 RGenAsp <--- FGenAsp
##    beta21  0.235458 0.119742   1.9664 4.9255e-02 FGenAsp <--- RGenAsp
##    lam21   1.062674 0.091967  11.5549 0.0000e+00 REdAsp <--- RGenAsp 
##    lam42   0.929727 0.071152  13.0668 0.0000e+00 FEdAsp <--- FGenAsp 
##    ps11    0.280987 0.046311   6.0674 1.2999e-09 RGenAsp <--> RGenAsp
##    ps22    0.263836 0.044902   5.8759 4.2067e-09 FGenAsp <--> FGenAsp
##    ps12   -0.022601 0.051649  -0.4376 6.6168e-01 FGenAsp <--> RGenAsp
##    theta1  0.412145 0.052211   7.8939 2.8866e-15 ROccAsp <--> ROccAsp
##    theta2  0.336148 0.053323   6.3040 2.9003e-10 REdAsp <--> REdAsp  
##    theta3  0.311194 0.046665   6.6687 2.5800e-11 FOccAsp <--> FOccAsp
##    theta4  0.404604 0.046733   8.6578 0.0000e+00 FEdAsp <--> FEdAsp  
##    
##     Iterations =  28 

# Fit the model using a numerical ram specification

ram.dhp <- matrix(c(
#               heads   to      from    param  start
                1,       1,     11,      0,     1,
                1,       2,     11,      1,     NA, # lam21
                1,       3,     12,      0,     1,
                1,       4,     12,      2,     NA, # lam42
                1,      11,      5,      3,     NA, # gam11
                1,      11,      6,      4,     NA, # gam12
                1,      11,      7,      5,     NA, # gam13
                1,      11,      8,      6,     NA, # gam14
                1,      12,      7,      7,     NA, # gam23
                1,      12,      8,      8,     NA, # gam24
                1,      12,      9,      9,     NA, # gam25
                1,      12,     10,     10,     NA, # gam26
                1,      11,     12,     11,     NA, # beta12
                1,      12,     11,     12,     NA, # beta21
                2,       1,      1,     13,     NA, # theta1
                2,       2,      2,     14,     NA, # theta2
                2,       3,      3,     15,     NA, # theta3
                2,       4,      4,     16,     NA, # theta4
                2,      11,     11,     17,     NA, # psi11
                2,      12,     12,     18,     NA, # psi22
                2,      11,     12,     19,     NA  # psi12
                ), ncol=5, byrow=TRUE)

params.dhp <- c('lam21', 'lam42', 'gam11', 'gam12', 'gam13', 'gam14',
                 'gam23',  'gam24',  'gam25',  'gam26',
                 'beta12', 'beta21', 'theta1', 'theta2', 'theta3', 'theta4',
                 'psi11', 'psi22', 'psi12')
                 
vars.dhp <- c('ROccAsp', 'REdAsp', 'FOccAsp', 'FEdAsp', 'RParAsp', 'RIQ',
                'RSES', 'FSES', 'FIQ', 'FParAsp', 'RGenAsp', 'FGenAsp')
                
sem.dhp.2 <- sem(ram.dhp, R.DHP, 329, params.dhp, vars.dhp, fixed.x=5:10)

summary(sem.dhp.2)
    
##     Model Chisquare =  26.697   Df =  15 Pr(>Chisq) = 0.031302
##     Chisquare (null model) =  872   Df =  45
##     Goodness-of-fit index =  0.98439
##     Adjusted goodness-of-fit index =  0.94275
##     RMSEA index =  0.048759   90 
##     Bentler-Bonnett NFI =  0.96938
##     Tucker-Lewis NNFI =  0.95757
##     Bentler CFI =  0.98586
##     BIC =  -60.244 
##    
##     Normalized Residuals
##       Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##    -0.7990 -0.1180  0.0000 -0.0120  0.0397  1.5700 
##    
##     Parameter Estimates
##           Estimate  Std Error z value Pr(>|z|)                       
##    lam21   1.062674 0.091967  11.5549 0.0000e+00 REdAsp <--- RGenAsp 
##    lam42   0.929727 0.071152  13.0668 0.0000e+00 FEdAsp <--- FGenAsp 
##    gam11   0.161224 0.038487   4.1890 2.8019e-05 RGenAsp <--- RParAsp
##    gam12   0.249653 0.044580   5.6001 2.1428e-08 RGenAsp <--- RIQ    
##    gam13   0.218404 0.043476   5.0235 5.0730e-07 RGenAsp <--- RSES   
##    gam14   0.071843 0.050335   1.4273 1.5350e-01 RGenAsp <--- FSES   
##    gam23   0.061894 0.051738   1.1963 2.3158e-01 FGenAsp <--- RSES   
##    gam24   0.228868 0.044495   5.1437 2.6938e-07 FGenAsp <--- FSES   
##    gam25   0.349039 0.044551   7.8346 4.6629e-15 FGenAsp <--- FIQ    
##    gam26   0.159535 0.040129   3.9755 7.0224e-05 FGenAsp <--- FParAsp
##    beta12  0.184226 0.096207   1.9149 5.5506e-02 RGenAsp <--- FGenAsp
##    beta21  0.235458 0.119742   1.9664 4.9255e-02 FGenAsp <--- RGenAsp
##    theta1  0.412145 0.052211   7.8939 2.8866e-15 ROccAsp <--> ROccAsp
##    theta2  0.336148 0.053323   6.3040 2.9003e-10 REdAsp <--> REdAsp  
##    theta3  0.311194 0.046665   6.6687 2.5800e-11 FOccAsp <--> FOccAsp
##    theta4  0.404604 0.046733   8.6578 0.0000e+00 FEdAsp <--> FEdAsp  
##    psi11   0.280987 0.046311   6.0674 1.2999e-09 RGenAsp <--> RGenAsp
##    psi22   0.263836 0.044902   5.8759 4.2067e-09 FGenAsp <--> FGenAsp
##    psi12  -0.022601 0.051649  -0.4376 6.6168e-01 RGenAsp <--> FGenAsp
    
     Iterations =  28 

# -------------------- Wheaton et al. alienation data ----------------------
    

S.wh <- matrix(c(
   11.834,     0,        0,        0,       0,        0,
    6.947,    9.364,     0,        0,       0,        0,
    6.819,    5.091,   12.532,     0,       0,        0,
    4.783,    5.028,    7.495,    9.986,    0,        0,
   -3.839,   -3.889,   -3.841,   -3.625,   9.610,     0,
  -21.899,  -18.831,  -21.748,  -18.775,  35.522,  450.288), 
  6, 6)
  
rownames(S.wh) <- colnames(S.wh) <-
    c('Anomia67','Powerless67','Anomia71','Powerless71','Education','SEI')
  
# This is the model in the SAS manual for PROC CALIS: A Recursive SEM with
# latent endogenous and exogenous variables.
# Curiously, both factor loadings for two of the latent variables are fixed.

model.wh.1 <- specify.model()
    Alienation67   ->  Anomia67,      NA,     1
    Alienation67   ->  Powerless67,   NA,     0.833
    Alienation71   ->  Anomia71,      NA,     1
    Alienation71   ->  Powerless71,   NA,     0.833 
    SES            ->  Education,     NA,     1     
    SES            ->  SEI,           lamb,   NA
    SES            ->  Alienation67,  gam1,   NA
    Alienation67   ->  Alienation71,  beta,   NA
    SES            ->  Alienation71,  gam2,   NA
    Anomia67       <-> Anomia67,      the1,   NA
    Anomia71       <-> Anomia71,      the1,   NA
    Powerless67    <-> Powerless67,   the2,   NA
    Powerless71    <-> Powerless71,   the2,   NA
    Education      <-> Education,     the3,   NA
    SEI            <-> SEI,           the4,   NA
    Anomia67       <-> Anomia71,      the5,   NA
    Powerless67    <-> Powerless71,   the5,   NA
    Alienation67   <-> Alienation67,  psi1,   NA
    Alienation71   <-> Alienation71,  psi2,   NA
    SES            <-> SES,           phi,    NA
    
                        
sem.wh.1 <- sem(model.wh.1, S.wh, 932)

summary(sem.wh.1)

##     Model Chisquare =  13.485   Df =  9 Pr(>Chisq) = 0.14186
##     Chisquare (null model) =  2131.4   Df =  15
##     Goodness-of-fit index =  0.99527
##     Adjusted goodness-of-fit index =  0.98896
##     RMSEA index =  0.023136   90 
##     Bentler-Bonnett NFI =  0.99367
##     Tucker-Lewis NNFI =  0.99647
##     Bentler CFI =  0.99788
##     BIC =  -48.051 
##    
##     Normalized Residuals
##        Min.  1st Qu.   Median     Mean  3rd Qu.     Max. 
##    -1.26000 -0.13100  0.00014 -0.02870  0.11400  0.87400 
##    
##     Parameter Estimates
##         Estimate  Std Error z value  Pr(>|z|)                                 
##    lamb   5.36880  0.433982  12.3710 0.0000e+00 SEI <--- SES                  
##    gam1  -0.62994  0.056128 -11.2233 0.0000e+00 Alienation67 <--- SES         
##    beta   0.59312  0.046820  12.6680 0.0000e+00 Alienation71 <--- Alienation67
##    gam2  -0.24086  0.055202  -4.3632 1.2817e-05 Alienation71 <--- SES         
##    the1   3.60787  0.200589  17.9864 0.0000e+00 Anomia67 <--> Anomia67        
##    the2   3.59494  0.165234  21.7567 0.0000e+00 Powerless67 <--> Powerless67  
##    the3   2.99366  0.498972   5.9996 1.9774e-09 Education <--> Education      
##    the4 259.57583 18.321121  14.1681 0.0000e+00 SEI <--> SEI                  
##    the5   0.90579  0.121710   7.4422 9.9032e-14 Anomia71 <--> Anomia67        
##    psi1   5.67050  0.422906  13.4084 0.0000e+00 Alienation67 <--> Alienation67
##    psi2   4.51481  0.334993  13.4773 0.0000e+00 Alienation71 <--> Alienation71
##    phi    6.61632  0.639506  10.3460 0.0000e+00 SES <--> SES                  
##    
##     Iterations =  78 

# The same model, but treating one loading for each latent variable as free.

model.wh.2 <- specify.model()
    Alienation67   ->  Anomia67,      NA,        1
    Alienation67   ->  Powerless67,   lamby,    NA
    Alienation71   ->  Anomia71,      NA,        1
    Alienation71   ->  Powerless71,   lamby,    NA 
    SES            ->  Education,     NA,        1     
    SES            ->  SEI,           lambx,    NA
    SES            ->  Alienation67,  gam1,     NA
    Alienation67   ->  Alienation71,  beta,     NA
    SES            ->  Alienation71,  gam2,     NA
    Anomia67       <-> Anomia67,      the1,     NA
    Anomia71       <-> Anomia71,      the1,     NA
    Powerless67    <-> Powerless67,   the2,     NA
    Powerless71    <-> Powerless71,   the2,     NA
    Education      <-> Education,     the3,     NA
    SEI            <-> SEI,           the4,     NA
    Anomia67       <-> Anomia71,      the5,     NA
    Powerless67    <-> Powerless71,   the5,     NA
    Alienation67   <-> Alienation67,  psi1,     NA
    Alienation71   <-> Alienation71,  psi2,     NA
    SES            <-> SES,           phi,      NA 

sem.wh.2 <- sem(model.wh.2, S.wh, 932)

summary(sem.wh.2)

##     Model Chisquare =  12.673   Df =  8 Pr(>Chisq) = 0.12360
##     Chisquare (null model) =  2131.4   Df =  15
##     Goodness-of-fit index =  0.99553
##     Adjusted goodness-of-fit index =  0.98828
##     RMSEA index =  0.025049   90 
##     Bentler-Bonnett NFI =  0.99405
##     Tucker-Lewis NNFI =  0.99586
##     Bentler CFI =  0.9978
##     BIC =  -42.026 
##    
##     Normalized Residuals
##         Min.   1st Qu.    Median      Mean   3rd Qu.      Max. 
##    -0.997000 -0.140000  0.000295 -0.028800  0.100000  0.759000 
##    
##     Parameter Estimates
##          Estimate  Std Error z value  Pr(>|z|)                                 
##    lamby   0.86261  0.033383  25.8402 0.0000e+00 Powerless67 <--- Alienation67 
##    lambx   5.35302  0.432591  12.3743 0.0000e+00 SEI <--- SES                  
##    gam1   -0.62129  0.056142 -11.0663 0.0000e+00 Alienation67 <--- SES         
##    beta    0.59428  0.047040  12.6335 0.0000e+00 Alienation71 <--- Alienation67
##    gam2   -0.23580  0.054684  -4.3121 1.6173e-05 Alienation71 <--- SES         
##    the1    3.74499  0.249823  14.9906 0.0000e+00 Anomia67 <--> Anomia67        
##    the2    3.49378  0.200754  17.4033 0.0000e+00 Powerless67 <--> Powerless67  
##    the3    2.97409  0.499661   5.9522 2.6456e-09 Education <--> Education      
##    the4  260.13252 18.298141  14.2163 0.0000e+00 SEI <--> SEI                  
##    the5    0.90377  0.121818   7.4190 1.1791e-13 Anomia71 <--> Anomia67        
##    psi1    5.47380  0.464073  11.7951 0.0000e+00 Alienation67 <--> Alienation67
##    psi2    4.36410  0.362722  12.0315 0.0000e+00 Alienation71 <--> Alienation71
##    phi     6.63576  0.640425  10.3615 0.0000e+00 SES <--> SES                  
##    
##     Iterations =  79 
##

# ----------------------- Thurstone data ---------------------------------------
#  Second-order confirmatory factor analysis, from the SAS manual for PROC CALIS

R.thur <- matrix(c(
        1.,       0,      0,      0,      0,      0,      0,      0,      0, 
         .828,   1.,      0,      0,      0,      0,      0,      0,      0, 
         .776,   .779,   1.,      0,      0,      0,      0,      0,      0, 
         .439,   .493,    .460,   1.,     0,      0,      0,      0,      0, 
         .432,   .464,    .425,   .674,   1.,     0,      0,      0,      0, 
         .447,   .489,    .443,   .590,    .541,  1.,     0,      0,      0, 
         .447,   .432,    .401,   .381,    .402,   .288,  1.,     0,      0, 
         .541,   .537,    .534,   .350,    .367,   .320,   .555,  1.,     0, 
         .380,   .358,    .359,   .424,    .446,   .325,   .598,   .452,  1.
            ), ncol=9, byrow=TRUE)
            
rownames(R.thur) <- colnames(R.Thur) <-
    c('Sentences','Vocabulary','Sent.Completion','First.Letters',
        '4.Letter.Words','Suffixes','Letter.Series','Pedigrees', 'Letter.Group')

model.thur <- specify.model()
    F1 -> Sentences,                      lam11, NA
    F1 -> Vocabulary,                     lam21, NA
    F1 -> Sent.Completion,                lam31, NA
    F2 -> First.Letters,                  lam41, NA
    F2 -> 4.Letter.Words,                 lam52, NA
    F2 -> Suffixes,                       lam62, NA
    F3 -> Letter.Series,                  lam73, NA
    F3 -> Pedigrees,                      lam83, NA
    F3 -> Letter.Group,                   lam93, NA
    F4 -> F1,                             gam1,  NA
    F4 -> F2,                             gam2,  NA
    F4 -> F3,                             gam3,  NA 
    Sentences <-> Sentences,              th1,   NA
    Vocabulary <-> Vocabulary,            th2,   NA
    Sent.Completion <-> Sent.Completion,  th3,   NA
    First.Letters <-> First.Letters,      th4,   NA
    4.Letter.Words <-> 4.Letter.Words,    th5,   NA
    Suffixes <-> Suffixes,                th6,   NA
    Letter.Series <-> Letter.Series,      th7,   NA
    Pedigrees <-> Pedigrees,              th8,   NA
    Letter.Group <-> Letter.Group,        th9,   NA
    F1 <-> F1,                            NA,     1
    F2 <-> F2,                            NA,     1
    F3 <-> F3,                            NA,     1
    F4 <-> F4,                            NA,     1

sem.thur <- sem(model.thur, R.thur, 213)

summary(sem.thur)

##     Model Chisquare =  38.196   Df =  24 Pr(>Chisq) = 0.033101
##     Chisquare (null model) =  1101.9   Df =  36
##     Goodness-of-fit index =  0.95957
##     Adjusted goodness-of-fit index =  0.9242
##     RMSEA index =  0.052822   90 
##     Bentler-Bonnett NFI =  0.96534
##     Tucker-Lewis NNFI =  0.98002
##     Bentler CFI =  0.98668
##     BIC =  -90.475 
##    
##     Normalized Residuals
##         Min.   1st Qu.    Median      Mean   3rd Qu.      Max. 
##    -9.72e-01 -4.16e-01 -4.20e-07  4.01e-02  9.39e-02  1.63e+00 
##    
##     Parameter Estimates
##          Estimate Std Error z value Pr(>|z|)                                       
##    lam11 0.51512  0.064964  7.9293  2.2204e-15 Sentences <--- F1                   
##    lam21 0.52031  0.065162  7.9849  1.3323e-15 Vocabulary <--- F1                  
##    lam31 0.48743  0.062422  7.8087  5.7732e-15 Sent.Completion <--- F1             
##    lam41 0.52112  0.063137  8.2538  2.2204e-16 First.Letters <--- F2               
##    lam52 0.49707  0.059673  8.3298  0.0000e+00 4.Letter.Words <--- F2              
##    lam62 0.43806  0.056479  7.7562  8.6597e-15 Suffixes <--- F2                    
##    lam73 0.45244  0.071371  6.3392  2.3100e-10 Letter.Series <--- F3               
##    lam83 0.41729  0.061037  6.8367  8.1020e-12 Pedigrees <--- F3                   
##    lam93 0.40763  0.064524  6.3175  2.6584e-10 Letter.Group <--- F3                
##    gam1  1.44381  0.264173  5.4654  4.6184e-08 F1 <--- F4                          
##    gam2  1.25383  0.216597  5.7888  7.0907e-09 F2 <--- F4                          
##    gam3  1.40655  0.279331  5.0354  4.7681e-07 F3 <--- F4                          
##    th1   0.18150  0.028400  6.3907  1.6517e-10 Sentences <--> Sentences            
##    th2   0.16493  0.027797  5.9334  2.9678e-09 Vocabulary <--> Vocabulary          
##    th3   0.26713  0.033468  7.9816  1.5543e-15 Sent.Completion <--> Sent.Completion
##    th4   0.30150  0.050686  5.9484  2.7073e-09 First.Letters <--> First.Letters    
##    th5   0.36450  0.052358  6.9617  3.3618e-12 4.Letter.Words <--> 4.Letter.Words  
##    th6   0.50642  0.059963  8.4455  0.0000e+00 Suffixes <--> Suffixes              
##    th7   0.39033  0.061599  6.3367  2.3474e-10 Letter.Series <--> Letter.Series    
##    th8   0.48137  0.065388  7.3618  1.8141e-13 Pedigrees <--> Pedigrees            
##    th9   0.50510  0.065227  7.7437  9.5479e-15 Letter.Group <--> Letter.Group      
##
##    Iterations =  54 
##

#------------------------- Kerchoff/Kenney path analysis ---------------------
# An observed-variable recursive SEM from the LISREL manual

R.kerch <- matrix(c(
    1,      0,      0,      0,      0,      0,      0,
    -.100,  1,      0,      0,      0,      0,      0,
     .277,  -.152,  1,      0,      0,      0,      0,
     .250,  -.108,  .611,   1,      0,      0,      0,
     .572,  -.105,  .294,   .248,   1,      0,      0,
     .489,  -.213,  .446,   .410,   .597,   1,      0,
     .335,  -.153,  .303,   .331,   .478,   .651,   1),
     ncol=7, byrow=TRUE)

rownames(R.kerch) <- colnames(R.kerch) <- c('Intelligence','Siblings',
    'FatherEd','FatherOcc','Grades','EducExp','OccupAsp')

    
model.kerch <- specify.model()
    Intelligence -> Grades,       gam51,    NA
    Siblings -> Grades,           gam52,    NA
    FatherEd -> Grades,           gam53,    NA
    FatherOcc -> Grades,          gam54,    NA
    Intelligence -> EducExp,      gam61,    NA
    Siblings -> EducExp,          gam62,    NA
    FatherEd -> EducExp,          gam63,    NA
    FatherOcc -> EducExp,         gam64,    NA
    Grades -> EducExp,            beta65,   NA
    Intelligence -> OccupAsp,     gam71,    NA
    Siblings -> OccupAsp,         gam72,    NA
    FatherEd -> OccupAsp,         gam73,    NA
    FatherOcc -> OccupAsp,        gam74,    NA
    Grades -> OccupAsp,           beta75,   NA
    EducExp -> OccupAsp,          beta76,   NA
    Grades <-> Grades,            psi5,     NA
    EducExp <-> EducExp,          psi6,     NA
    OccupAsp <-> OccupAsp,        psi7,     NA

                       
sem.kerch <- sem(model.kerch, R.kerch, 737, fixed.x=c('Intelligence','Siblings',
    'FatherEd','FatherOcc'))
    
summary(sem.kerch)

##     Model Chisquare =  3.2685e-13   Df =  0 Pr(>Chisq) = NA
##     Chisquare (null model) =  1664.3   Df =  21
##     Goodness-of-fit index =  1
##     BIC =  3.2685e-13 
##    
##     Normalized Residuals
##         Min.   1st Qu.    Median      Mean   3rd Qu.      Max. 
##    -4.28e-15  0.00e+00  0.00e+00  7.62e-16  1.47e-15  5.17e-15 
##    
##     Parameter Estimates
##           Estimate  Std Error z value  Pr(>|z|)                             
##    gam51   0.525902 0.031182  16.86530 0.0000e+00 Grades <--- Intelligence  
##    gam52  -0.029942 0.030149  -0.99314 3.2064e-01 Grades <--- Siblings      
##    gam53   0.118966 0.038259   3.10951 1.8740e-03 Grades <--- FatherEd      
##    gam54   0.040603 0.037785   1.07456 2.8257e-01 Grades <--- FatherOcc     
##    gam61   0.160270 0.032710   4.89979 9.5940e-07 EducExp <--- Intelligence 
##    gam62  -0.111779 0.026876  -4.15899 3.1966e-05 EducExp <--- Siblings     
##    gam63   0.172719 0.034306   5.03461 4.7882e-07 EducExp <--- FatherEd     
##    gam64   0.151852 0.033688   4.50758 6.5571e-06 EducExp <--- FatherOcc    
##    beta65  0.405150 0.032838  12.33799 0.0000e+00 EducExp <--- Grades       
##    gam71  -0.039405 0.034500  -1.14215 2.5339e-01 OccupAsp <--- Intelligence
##    gam72  -0.018825 0.028222  -0.66700 5.0477e-01 OccupAsp <--- Siblings    
##    gam73  -0.041333 0.036216  -1.14126 2.5376e-01 OccupAsp <--- FatherEd    
##    gam74   0.099577 0.035446   2.80924 4.9658e-03 OccupAsp <--- FatherOcc   
##    beta75  0.157912 0.037443   4.21738 2.4716e-05 OccupAsp <--- Grades      
##    beta76  0.549593 0.038260  14.36486 0.0000e+00 OccupAsp <--- EducExp     
##    psi5    0.650995 0.033946  19.17743 0.0000e+00 Grades <--> Grades        
##    psi6    0.516652 0.026943  19.17590 0.0000e+00 EducExp <--> EducExp      
##    psi7    0.556617 0.029026  19.17644 0.0000e+00 OccupAsp <--> OccupAsp    
##    
##     Iterations =  0 

#------------------- McArdle/Epstein latent-growth-curve model -----------------
# This model, from McArdle and Epstein (1987, p.118), illustrates the use of a 
# raw moment matrix to fit a model with an intercept. (The example was suggested
# by Mike Stoolmiller.)

M.McArdle <- scan()
365.661     0           0           0           0
503.175     719.905     0           0           0
675.656     958.479    1303.392     0           0
890.680    1265.846    1712.475    2278.257     0
 18.034      25.819      35.255      46.593     1.000
 
M.McArdle <- matrix(M.McArdle, 5, 5, byrow=TRUE)
rownames(M.McArdle) <- colnames(M.McArdle) <- scan(what="")
WISC1 WISC2 WISC3 WISC4 UNIT
 
mod.McArdle <- specify.model()
C -> WISC1, NA, 6.07
C -> WISC2, B2, NA
C -> WISC3, B3, NA
C -> WISC4, B4, NA
UNIT -> C, Mc, NA
C <-> C, Vc, NA,
WISC1 <-> WISC1, Vd, NA
WISC2 <-> WISC2, Vd, NA
WISC3 <-> WISC3, Vd, NA
WISC4 <-> WISC4, Vd, NA

sem.McArdle <- sem(mod.McArdle, M.McArdle, 204, fixed.x="UNIT", raw=TRUE)
summary(sem.McArdle)

##    Model fit to raw moment matrix.
##    
##     Model Chisquare =  83.791   Df =  8 Pr(>Chisq) = 8.4377e-15
##     BIC =  41.246 
##    
##     Normalized Residuals
##        Min.  1st Qu.   Median     Mean  3rd Qu.     Max. 
##    -0.15300 -0.01840  0.00132 -0.00576  0.02400  0.07760 
##    
##     Parameter Estimates
##       Estimate Std Error z value Pr(>|z|)                 
##    B2  8.61354 0.135438  63.5976 0        WISC2 <--- C    
##    B3 11.64054 0.168854  68.9387 0        WISC3 <--- C    
##    B4 15.40323 0.213071  72.2916 0        WISC4 <--- C    
##    Mc  3.01763 0.060690  49.7219 0        C <--- UNIT     
##    Vc  0.44343 0.047704   9.2955 0        C <--> C        
##    Vd 11.78832 0.674060  17.4885 0        WISC1 <--> WISC1
##    
##     Iterations =  37 

    ## End(Not run)

Besides, you may draw your path diagram (see the function path.diagram()) via the graph-drawing program dot; see Koutsofios and North (2002) and http://www.graphviz.org/ (also free & open-source). Below is a diagram I draw in dot:

In the end, I believe the function _boot.sem() _is worth our attention. As the model is too complicated, it’s natural to doubt its stability, so _bootstrapping _surely will be a good choice, which might be more reliable than the coefficients’ standard errors directly computed under the assumption of multi-normality. Bear in mind that bootstrapping takes some time, so be patient.