sem                   package:sem                   R Documentation

_G_e_n_e_r_a_l _S_t_r_u_c_t_u_r_a_l _E_q_u_a_t_i_o_n _M_o_d_e_l_s

_D_e_s_c_r_i_p_t_i_o_n:

     '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.

_U_s_a_g_e:

     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, ...)

_A_r_g_u_m_e_n_t_s:

     ram: RAM specification, which is a simple encoding of the path
          diagram for the model. The 'ram' matrix may be given either
          in symbolic  form (as a 'mod' object, as returned by the
          'specify.model' function, or as a character matrix), invoking
          'sem.mod', which calls 'sem.default' after setting up the
          model, or (less conveniently) in numeric form, invoking
          'sem.default' directly (see *Details* below).

       S: covariance matrix among observed variables; may be input as a
          symmetric matrix, or as a lower- or upper-triangular matrix.
          'S' may also be a raw (i.e., ``uncorrected'') moment matrix -
          that is, a sum-of-squares-and-products matrix divided by 'N'.
          This form of input is useful for fitting models with
          intercepts, in which case the moment matrix should include
          the mean square and cross-products for a unit variable all of
          whose entries are 1; of course, the raw mean square for the
          unit variable is 1. Raw-moment matrices may be computed by
          'raw.moments'.

       N: number of observations on which the covariance matrix is
          based.

obs.variables: names of observed variables, by default taken from the
          row names of the covariance matrix 'S'.

 fixed.x: names (if the 'ram' matrix is given in symbolic form) or
          indices  (if it is in numeric form) of fixed exogenous
          variables. Specifying these obviates the necessity of having
          to fix the variances and covariances among these  variables
          (and produces correct degrees of freedom for the model
          chisquare).

     raw: 'TRUE' if 'S' is a raw moment matrix, as opposed to a
          covariance matrix; the default is 'FALSE'.

   debug: if 'TRUE', some information is printed to help you debug the
          symbolic model specification; for example, if a variable name
          is misspelled, 'sem' will assume that the variable is a (new)
          latent variable. The default is 'FALSE'.

     ...: arguments to be passed down, including from 'sem.default' to
          the 'nlm' optimizer.

param.names: names of the t free parameters, given in their numerical
          order; default names are 'Param1', ..., 'Param'_t_. Note:
          Should not be specified when the 'ram' matrix is given in
          symbolic form.

var.names: names of the m entries of the v vector (typically the
          observed and latent variables - see below), given in their
          numerical order; default names are 'Var1', ..., 'Var'_m_.
          Note: Should not be specified when the 'ram' matrix is given
          in symbolic form.

analytic.gradient: if 'TRUE' (the default), then analytic first
          derivatives are used in the maximization of the likelihood;
          otherwise numeric derivatives are used.

    warn: if 'TRUE', warnings produced by the optimization function
          will be printed. This should generally not be necessary,
          since 'sem' prints its own warning, and saves information
          about convergence. The default is 'FALSE'.

 maxiter: the maximum number of iterations for the optimization
          performed by the 'nlm' function, to be passed to it via its
          'iterlim' argument.

par.size: the anticipated size of the free parameters; if '"ones"', a
          vector of ones is used; if '"startvalues"', taken from the
          start values. You can try changing this argument if you
          encounter convergence problems. The default is
          '"startvalues"' if the largest input variance is at least 100
          times the smallest, and '"ones"' otherwise.

   refit: if 'TRUE' (the default), attempt to refit the model
          eliminating apparently aliased parameters if 
          under-identification is detected.

start.tol, tol: if the magnitude of an automatic start value is less
          than 'start.tol', then it is set to 'start.tol'; defaults to
          1E-6.

object, x: an object of class 'sem' returned by the 'sem' function.

  digits: number of digits for printed output.

conf.level: level for confidence interval for the RMSEA index (default
          is .9).

_D_e_t_a_i_l_s:

     The model is set up using RAM (`reticular action model' - don't
     ask!)  notation - a simple format for specifying general
     structural equation models by coding the `arrows' in the path
     diagram for the model (see, e.g., McArdle and McDonald, 1984).

     The variables in the v vector in the model (typically, the
     observed and unobserved variables, but not error variables) are
     numbered from 1 to m.   the RAM matrix contains one row for each
     (free or constrained) parameter of the model, and may be specified
     either in symbolic format or in numeric format.

     A symbolic 'ram' matrix consists of three columns, as follows:

     _1. _A_r_r_o_w _s_p_e_c_i_f_i_c_a_t_i_o_n: This is a simple formula, of the form '"A
          -> B"' or, equivalently, '"B <- A"' for a regression
          coefficient (i.e., a single-headed or directional arrow); 
          '"A <-> A"' for a variance or '"A <-> B"' for a covariance
          (i.e., a double-headed or bidirectional arrow). Here, 'A' and
          'B' are variable names in the model. If a name does not
          correspond to an observed variable, then it is assumed to be
          a latent variable. Spaces can appear freely in an arrow
          specification, and there can be any number of hyphens in the
          arrows, including zero: Thus, e.g., '"A->B"', '"A --> B"',
          and '"A>B"' are all legitimate and equivalent.

     _2. _P_a_r_a_m_e_t_e_r _n_a_m_e: The name of the regression coefficient,
          variance, or covariance specified by the arrow. Assigning the
          same name to two or more arrows results in an equality
          constraint. Specifying the parameter name as 'NA' produces a
          fixed parameter.

     _3. _V_a_l_u_e: start value for a free parameter or value of a fixed
          parameter. If given as 'NA', 'sem' will compute the start
          value.

     It is simplest to construct the RAM matrix with the
     'specify.model' function, which returns an object of class 'mod'.
     This process is illustrated in the examples below. 

     A numeric 'ram' matrix consists of five columns, as follows:

     _1. _N_u_m_b_e_r _o_f _a_r_r_o_w _h_e_a_d_s: 1 (directed arrow) or 2 (covariance).

     _2. _A_r_r_o_w __t_o_: index of the variable at the head of a directional
          arrow, or at one end of a bidirectional arrow. Observed
          variables should be assigned the numbers 1 to n, where n is
          the number of rows/columns in the covariance matrix 'S', with
          the indices corresponding to the variables' positions in 'S'.
          Variable indices above n represent latent variables.

     _3. _A_r_r_o_w __f_r_o_m_: the index of the variable at the tail of  a
          directional arrow, or at the other end of a bidirectional
          arrow.

     _4. _P_a_r_a_m_e_t_e_r _n_u_m_b_e_r: free parameters are numbered from 1 to t, but
          do not necessarily appear in consecutive order. Fixed
          parameters are given the number 0. Equality contraints are
          specified by assigning two or more  parameters the same
          number.

     _5. _V_a_l_u_e: start value for a free parameter, or value of a fixed
          parameter. If given as 'NA', the program will compute a start
          value, by a slight modification of the method described by
          McDonald and Hartmann (1992). _Note:_ In some circumstances,
          some start values are selected randomly; this might produce
          small differences in the parameter estimates when the program
          is rerun.

     'sem' fits the model by calling the 'nlm' optimizer to minimize
     the negative log-likelihood for the model.  If 'nlm' fails to
     converge, a warning message is printed.

     The RAM formulation of the general structural equation model is
     given by the basic equation

                              v = Av + u

     where v and u are vectors of random variables (observed or
     unobserved), and the parameter matrix A contains regression
     coefficients, symbolized by single-headed arrows in a path
     diagram. Another parameter matrix, 

                              P = E(uu')

     contains covariances among the elements of u (assuming that the
     elements of u have zero means). Usually v contains endogenous and
     exogenous observed and unobserved variables, but not error
     variables (see the examples below).

     The 'startvalues' function may be called directly, but is usually
     called by 'sem.default'.

_V_a_l_u_e:

     'sem' returns an object of class 'sem', with the following
     elements: 

     ram: RAM matrix, including any rows generated for covariances
          among fixed exogenous variables; column 5 includes computed
          start values.

   coeff: estimates of free parameters.

criterion: fitting criterion - negative log-liklihood divided by N - 1
          (or N if the model is fit to a raw moment matrix).

     cov: estimated asymptotic covariance matrix of parameter
          estimates.

       S: observed covariance matrix.

       J: RAM selection matrix, J, which picks out observed variables.

       C: model-reproduced covariance matrix.

       A: RAM A matrix.

       P: RAM P matrix.

   n.fix: number of fixed exogenous variables.

       n: number of observed variables.

       N: number of observations.

       m: number of variables (observed plus unobserved).

       t: number of free parameters.

par.posn: indices of free parameters.

var.names: vector of variable names.

observed: indices of observed variables.

convergence: convergence code returned by 'nlm' (a code > 2 indicates a
          problem).

iterations: number of iterations performed.

     raw: 'TRUE' if the model is fit to a raw moment matrix, 'FALSE'
          otherwise.

_W_a_r_n_i_n_g:

     A common error is to fail to specify variance or covariance terms
     in the model, which are denoted  by double-headed arrows, '<->'.

     In general, every observed or latent variable in the model should
     be associated  with a variance or error variance.  This may be a
     free parameter to estimate or a fixed constant (as in the case of
     a latent exogenous  variable for which you wish to fix the
     variance, e.g., to 1). Again in general, there will be an _error
     variance_ associated with each  endogenous variable in the model
     (i.e., each variable to which at least one single-headed arrow
     points - including observed indicators of latent variables),  and
     a _variance_ associated with each exogenous variable (i.e., each
     variable that appears only at the tail of single-headed arrows,
     never at the head).

     To my knowledge, the only _apparent_ exception to this rule is for
     observed variables that  are declared to be fixed exogenous
     variables. In this case, the program generates the necessary
     (fixed-constant) variances and covariances automatically.

     If there are missing variances, a warning message will be printed,
     and estimation will almost surely  fail in some manner. Missing
     variances might well indicate that there are missing covariances
     too, but it is not possible  to deduce this in a mechanical
     manner.

_A_u_t_h_o_r(_s):

     John Fox jfox@mcmaster.ca

_R_e_f_e_r_e_n_c_e_s:

     Bollen, K. A. (1989)  _Structural Equations With Latent
     Variables._ Wiley.

     McArdle, J. J. and Epstein, D. (1987) Latent growth curves within
     developmental structural equation models. _Child Development_
     *58*, 110-133.

     McArdle, J. J. and McDonald, R. P. (1984) Some algebraic
     properties of the reticular action model. _British Journal of
     Mathematical and Statistical Psychology_ *37*, 234-251.

     McDonald, R. P. and Hartmann, W. M. (1992) A procedure for
     obtaining initial values of parameters in the RAM model.
     _Multivariate Behavioral Research_ *27*, 57-76.

     Raftery, A. E. (1993) Bayesian model selection in structural
     equation models. In Bollen, K. A. and Long, J. S. (eds.) _Testing
     Structural Equation Models_, Sage.

_S_e_e _A_l_s_o:

     'raw.moments', 'nlm'

_E_x_a_m_p_l_e_s:

     # ------------- 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)

     # 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

     obs.vars.dhp <- c('ROccAsp', 'REdAsp', 'FOccAsp', 'FEdAsp', 'RParAsp', 'RIQ',
                     'RSES', 'FSES', 'FIQ', 'FParAsp')

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

     summary(sem.dhp.1)

     ##  Model Chisquare =  26.697   Df =  15 Pr(>Chisq) = 0.031302
     ##  Goodness-of-fit index =  0.98439
     ##  Adjusted goodness-of-fit index =  0.94275
     ##  RMSEA index =  0.048759   90 
     ##  BIC =  -94.782 
     ## 
     ##  Normalized Residuals
     ##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
     ## -0.8010 -0.1180  0.0000 -0.0120  0.0398  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.2068e-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 

     ## Not run: 

     # 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
     ##  Goodness-of-fit index =  0.98439
     ##  Adjusted goodness-of-fit index =  0.94275
     ##  RMSEA index =  0.048759   90 
     ##  BIC =  -94.782 
     ## 
     ##  Normalized Residuals
     ##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
     ## -0.8010 -0.1180  0.0000 -0.0120  0.0398  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)
       
     # 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
         
                             
     obs.vars.wh <- c('Anomia67','Powerless67','Anomia71','Powerless71','Education','SEI')

     sem.wh.1 <- sem(model.wh.1, S.wh, 932, obs.vars.wh)

     summary(sem.wh.1)

     ##  Model Chisquare =  13.485   Df =  9 Pr(>Chisq) = 0.14186
     ##  Goodness-of-fit index =  0.99527
     ##  Adjusted goodness-of-fit index =  0.98896
     ##  RMSEA index =  0.023136   90 
     ##  BIC =  -64.177 
     ## 
     ##  Normalized Residuals
     ##     Min.  1st Qu.   Median     Mean  3rd Qu.     Max. 
     ## -1.26000 -0.13100  0.00014 -0.02870  0.11400  0.87500 
     ## 
     ##  Parameter Estimates
     ##      Estimate  Std Error z value  Pr(>|z|)                                 
     ## lamb   5.36880  0.433981  12.3710 0.0000e+00 SEI <--- SES                  
     ## gam1  -0.62994  0.056128 -11.2234 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.498971   5.9997 1.9772e-09 Education <--> Education      
     ## the4 259.57583 18.321099  14.1681 0.0000e+00 SEI <--> SEI                  
     ## the5   0.90579  0.121711   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, obs.vars.wh)

     summary(sem.wh.2)

     ##  Model Chisquare =  12.673   Df =  8 Pr(>Chisq) = 0.12360
     ##  Goodness-of-fit index =  0.99553
     ##  Adjusted goodness-of-fit index =  0.98828
     ##  RMSEA index =  0.025049   90 
     ##  BIC =  -56.36 
     ## 
     ##  Normalized Residuals
     ##      Min.   1st Qu.    Median      Mean   3rd Qu.      Max. 
     ## -0.998000 -0.140000  0.000296 -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.499662   5.9522 2.6457e-09 Education <--> Education      
     ## the4  260.13252 18.298143  14.2163 0.0000e+00 SEI <--> SEI                  
     ## the5    0.90377  0.121818   7.4190 1.1791e-13 Anomia71 <--> Anomia67        
     ## psi1    5.47380  0.464074  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)

     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

                             
     obs.vars.thur <- c('Sentences','Vocabulary','Sent.Completion','First.Letters',
                     '4.Letter.Words','Suffixes','Letter.Series','Pedigrees',
                     'Letter.Group')

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

     summary(sem.thur)

     ##  Model Chisquare =  38.196   Df =  24 Pr(>Chisq) = 0.033101
     ##  Goodness-of-fit index =  0.95957
     ##  Adjusted goodness-of-fit index =  0.9242
     ##  RMSEA index =  0.052822   90 
     ##  BIC =  -143.21 
     ## 
     ##  Normalized Residuals
     ##      Min.   1st Qu.    Median      Mean   3rd Qu.      Max. 
     ## -9.75e-01 -4.17e-01 -2.90e-06  4.02e-02  9.41e-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.071372  6.3392  2.3103e-10 Letter.Series <--- F3               
     ## lam83 0.41729  0.061037  6.8367  8.1031e-12 Pedigrees <--- F3                   
     ## lam93 0.40763  0.064524  6.3175  2.6587e-10 Letter.Group <--- F3                
     ## gam1  1.44381  0.264173  5.4654  4.6186e-08 F1 <--- F4                          
     ## gam2  1.25383  0.216597  5.7888  7.0904e-09 F2 <--- F4                          
     ## gam3  1.40655  0.279332  5.0354  4.7684e-07 F3 <--- F4                          
     ## th1   0.18150  0.028401  6.3907  1.6517e-10 Sentences <--> Sentences            
     ## th2   0.16493  0.027797  5.9334  2.9679e-09 Vocabulary <--> Vocabulary          
     ## th3   0.26713  0.033469  7.9816  1.5543e-15 Sent.Completion <--> Sent.Completion
     ## th4   0.30150  0.050686  5.9484  2.7074e-09 First.Letters <--> First.Letters    
     ## th5   0.36450  0.052358  6.9617  3.3618e-12 4.Letter.Words <--> 4.Letter.Words  
     ## th6   0.50641  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 =  53 

     #------------------------- 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
     ##  Goodness-of-fit index =  1
     ##  Adjusted goodness-of-fit index =  NA
     ##  RMSEA index =  Inf   90 
     ##  BIC =  NA 
     ## 
     ##  Normalized Residuals
     ##      Min.   1st Qu.    Median      Mean   3rd Qu.      Max. 
     ## -4.26e-15 -1.35e-15  0.00e+00 -4.17e-16  0.00e+00  1.49e-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.7881e-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
     ##  Goodness-of-fit index =  0.84537
     ##  Adjusted goodness-of-fit index =  0.71007
     ##  RMSEA index =  0.2155   90 
     ##  BIC =  28.371 
     ## 
     ##  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.674059  17.4885 0        WISC1 <--> WISC1
     ##
     ##  Iterations =  37 
         ## End(Not run)

