Factor analysis with -factor- and -sem-

use http://www.stata-press.com/data/r15/bg2, clear

factor bg2cost1 bg2cost5 bg2cost6, ml factors(1)
sem (Y -> bg2cost1 bg2cost5 bg2cost6), standardized

Mediation analysis

This replicates two illustrations from Iacobucci (2008).

// Section 3.3 (pp. 21-23)
  // Illustration in section 3.4.1 (pp. 25-26) is identical, just the number of 
  // observations should be 100

  // Read in data
clear
ssd init y x m
ssd set observations 50
#delimit ;
ssd set corr 
 1.00 \
 0.45 1.00 \
 0.63 0.55 1.00
;
#delimit cr

// SEM approach
sem (m <- x) (y <- m x)
  // Calculation by hand
di (_b[m:x] * _b[y:m]) / ((_b[m:x] * _b[y:m]) +  _b[y:x])
  // Calculation via estat, teffects
estat teffects, compact
matrix b_indirect = r(indirect)
matrix b_total = r(total)

scalar indirect  = el(b_indirect, 1, 3)
scalar total = el(b_total, 1, 3)

di "Proportion of total effect mediated by M: "  indirect/total
  
// Section 4.2 (pp. 35-38)
  // Read in data
clear
ssd init y x m q
ssd set observations 50
#delimit ;
ssd set corr 
 1.00 \
 0.45 1.00 \
 0.63 0.55 1.00 \
 0.40 0.40 0.40 1.00
;
#delimit cr

// Figure 4.3
  // Baseline
sem (m <- x) (y <- m x) 
  // (a)
sem (m <- x) (y <- m x) (x <- q)
  // (b)
sem (m <- x) (y <- m x) (x -> q)
  // (d)
sem (m <- x) (y <- m x) (m -> q)
  // (e)
sem (m <- x) (y <- m x) (y -> q)

// Figure 4.4
  // (c) r's = .40
sem (m <- x) (y <- m x) (m <- q)  
  // (e) r's = .40
sem (m <- x) (y <- m x) (y <- q)  

clear
ssd init y x m q
ssd set observations 50
#delimit ;
ssd set corr 
 1.00 \
 0.45 1.00 \
 0.63 0.55 1.00 \
 0.70 0.70 0.70 1.00
;
#delimit cr
  // (c') r's = .70
sem (m <- x) (y <- m x) (m <- q)  
  // (e') r's = .70
sem (m <- x) (y <- m x) (y <- q)  

//

Reference

Iacobucci, Dawn. 2008. Mediation Analysis. Sage. doi: 10.4135/9781412984966

Analysis of incomplete data with full information ML using Stata

The code below allows replicating the example of Allison (2002, pp. 25-27).
  

// Table 4.6
use "https://statisticalhorizons.com/wp-content/uploads/college.dta", clear



eststo clear
eststo: sem (gradrat act <- csat lenroll private stufac rmbrd), cov(e.gradrat*e.act) method(mlmv)

#delimit ;
esttab using test.tex, cells("b(fmt(3) label(Coefficient))
                              se(fmt(3) label(Standard Error)) 
                              t(fmt(2) label(t Statistic)) 
                              p(fmt(4) label(p Value))")
                       order(_cons) coeflabel(_cons "Intercept")
                       nomtitle nonumber title(Regression that predicts GRADRAT Using Direct ML)
                       keep(gradrat:) eqlabels("", none) // Removes equation label
                       booktabs replace;
#delimit cr

Reference

Allison, Paul D. 2002. Missing Data. Sage. doi: 10.4135/9781412985079

IV regression using the -sem- command

// Read in data from Angrist and Krueger (1991)
// https://economics.mit.edu/faculty/angrist/data1/data/angkru1991
infile lwklywge educ yob qob pob using asciiqob.txt, clear

// Generate dummy variables as SEM command does not take factor variables
qui tabulate qob, gen(qobx)
qui tabulate yob, gen(yobx)
qui tabulate pob, gen(pobx)
drop qobx1 yobx1 pobx1      // Get rid of reference category

eststo clear

// Model 2 of Table 4.1.1 of Mostly Harmless Econometrics
eststo: regress lwklywge educ yobx* pobx*, robust                              

// Model 6 of Table 4.1.1 of Mostly Harmless Econometrics
eststo: ivregress 2sls lwklywge yobx* pobx* (educ = qobx*), robust

// Model 6 of Table 4.1.1 using the sem command
eststo: sem (lwklywge <- yobx* pobx* educ) (educ <- pobx* qobx*), cov(e.lwklywge*e.educ) 


esttab, b(3) se(3) nostar drop(_cons educ:) ///
        indicate("9 year-of-birth dummies = yobx*" ///
                 "50 state-of-birth dummies = pobx*") ///
        title("OLS and 2SLS estimates of the economic returns to schooling") ///
        coeflabel(educ "Years of education") ///
        mtitles("OLS" "-ivregress-" "-sem-") nonumbers varwidth(30) ///
        eqlabels("", none) // Removes equation label

Chapter 6 of Allison's (2009) book on Fixed-Effects Regression Models using Stata

use http://statisticalhorizons.com/wp-content/uploads/nlsy.dta, clear
 
// Random effects model
sem (Falpha -> anti90@1 anti92@1 anti94@1) /// 
    (anti90 <- pov90@a self90@b black@c hispanic@d childage@e married@f gender@g momage@h momwork@i) /// 
    (anti92 <- pov92@a self92@b black@c hispanic@d childage@e married@f gender@g momage@h momwork@i) ///  
    (anti94 <- pov94@a self94@b black@c hispanic@d childage@e married@f gender@g momage@h momwork@i) /// 
   , cov(Falpha*pov90@0   Falpha*pov92@0    Falpha*pov94@0    /// 
         Falpha*self90@0  Falpha*self92@0   Falpha*self94@0   /// 
         Falpha*black@0   Falpha*hispanic@0 Falpha*childage@0 /// 
         Falpha*married@0 Falpha*gender@0   Falpha*momage@0   /// 
         Falpha*momwork@0)                                    /// 
     var(e.anti90@j       e.anti92@j        e.anti94@j)
estimates store random
  
// Fixed effects model
sem (Falpha -> anti90@1 anti92@1 anti94@1) /// 
    (anti90 <- pov90@a self90@b black@c hispanic@d childage@e married@f gender@g momage@h momwork@i) /// 
    (anti92 <- pov92@a self92@b black@c hispanic@d childage@e married@f gender@g momage@h momwork@i) ///  
    (anti94 <- pov94@a self94@b black@c hispanic@d childage@e married@f gender@g momage@h momwork@i) /// 
   , cov(Falpha*black@0   Falpha*hispanic@0 Falpha*childage@0 /// 
         Falpha*married@0 Falpha*gender@0   Falpha*momage@0   /// 
         Falpha*momwork@0)                                    /// 
     var(e.anti90@j       e.anti92@j        e.anti94@j)
estimates store fixed
  
// Compromise model
sem (Falpha -> anti90@1 anti92@1 anti94@1) /// 
    (anti90 <- pov90@a self90@b black@c hispanic@d childage@e married@f gender@g momage@h momwork@i) /// 
    (anti92 <- pov92@a self92@b black@c hispanic@d childage@e married@f gender@g momage@h momwork@i) ///  
    (anti94 <- pov94@a self94@b black@c hispanic@d childage@e married@f gender@g momage@h momwork@i) /// 
   , cov(Falpha*self90@0  Falpha*self92@0   Falpha*self94@0   /// 
         Falpha*black@0   Falpha*hispanic@0 Falpha*childage@0 /// 
         Falpha*married@0 Falpha*gender@0   Falpha*momage@0   /// 
         Falpha*momwork@0)                                    /// 
     var(e.anti90@j       e.anti92@j        e.anti94@j)
estimates store compromise

// Table 6.1
estimates table random fixed compromise, ///
          b(%9.3f) se(%9.3f) ///
          keep(anti90:self90   anti90:pov90    anti90:black ///
               anti90:hispanic anti90:childage anti90:married ///
               anti90:gender   anti90:momage   anti90:momwork)

// Table 6.2
// Shows covariances instead of correlations
estimates table fixed, b t keep(cov(pov90,Falpha):_cons cov(pov92,Falpha):_cons cov(pov94,Falpha):_cons ///
                            cov(self90,Falpha):_cons cov(self92,Falpha):_cons cov(self94,Falpha):_cons)

// Table 6.3
use http://statisticalhorizons.com/wp-content/uploads/occ.dta, clear  
// Thanks to http://www.stata.com/statalist/archive/2013-10/msg00019.html

sem (Falpha -> pf2@1 pf3@1 pf4@1) ///
    (pf4 <- pf3@a mdwgf3@b) /// 
    (pf3 <- pf2@a mdwgf2@b) /// 
    (pf2 <- pf1@a mdwgf1@b ERR1@1), ///
    cov(e.pf2@0) cov(ERR1*_oexogenous@0 ERR1*Falpha@0 ERR1*mdwgf3) method(mlmv)
estimates store medianwage

sem (Falpha -> mdwgf2@1 mdwgf3@1 mdwgf4@1) /// 
    (mdwgf4 <- pf3@a mdwgf3@b) /// 
    (mdwgf3 <- pf2@a mdwgf2@b) /// 
    (mdwgf2 <- pf1@a mdwgf1@b ERR1@1), ///
    cov(e.mdwgf2@0) cov(ERR1*_oexogenous@0 ERR1*Falpha@0 ERR1*pf3) method(mlmv)
estimates store proportionfemale

esttab medianwage proportionfemale, ///
       b(%9.3f) se(%9.3f) ///
       keep(main:pf1 main:mdwgf1)

Reference

Allison, Paul D. 2009. Fixed Effects Regression Models. Sage. doi: 10.4135/9781412993869

Preacher et al. (2008): Latent Growth Curve Modeling in Stata

The website accompanying Preacher et al.'s (2008) book on latent growth curve modeling provide syntax files for Lisrel, Mplus, and Mx, but not for Stata. These are the Stata commands for replicating the first four models presented in the book.

version 13

// Chapter 2

clear
// Read in data from clsn_cov.dat
ssd init clsn1 clsn3 clsn4 clsn5 clsn6 sex
ssd set observations 851
ssd set means 37.9542 37.2785 37.0463 36.5696 36.1363 0.49
#delimit ;
ssd set cov 
 6.3944 \
 3.2716 7.5282 \
 4.1435 6.0804 10.7290 \
 3.7058 5.1597 6.5672 10.2920 \
 4.1286 5.7608 7.2365 7.6463 12.9085 \
-0.0940 -0.0390 -0.1521 -0.1104 -0.1469 0.2502
;
#delimit cr

// Table 2.2
ssd list

// Model 0: The null model
// Only mean intercept and residual variance are estimated

sem (Intercept@1 -> clsn1, ) ///  // Set loadings to 1 
    (Intercept@1 -> clsn3, ) ///  // for estimating constant
    (Intercept@1 -> clsn4, ) ///
    (Intercept@1 -> clsn5, ) ///
    (Intercept@1 -> clsn6, ) ///
    (clsn1 <- _cons@a, ) ///      // Constrain means to be
    (clsn3 <- _cons@a, ) ///      // the same
    (clsn4 <- _cons@a, ) ///
    (clsn5 <- _cons@a, ) ///
    (clsn6 <- _cons@a, ) ///
  , latent(Intercept ) ///
    cov(Intercept@0 ///  Constrain intercept variance to 0
        e.clsn1@b e.clsn3@b ///  // Constrain residual variances
        e.clsn4@b e.clsn5@b ///  // to be the same
        e.clsn6@b) ///
    nocapslatent
estimates store m0
estat gof, stats(chi2 rmsea indices residuals)
// Non-normed fit index (NNFI) is called 
// Tucker-Lewis index (TLI) in Stata


// Model 1: Random intercept model (Table 2.3)
// Only mean intercept, intercept variance, and residual variance are 
// estimated



sem (Intercept@1 -> clsn1, ) ///  // Set loadings to 1 
    (Intercept@1 -> clsn3, ) ///  // for estimating constant
 (Intercept@1 -> clsn4, ) ///
 (Intercept@1 -> clsn5, ) ///
 (Intercept@1 -> clsn6, ) ///
 (clsn1 <- _cons@a, ) ///      // Constrain means to be
 (clsn3 <- _cons@a, ) ///      // the same
 (clsn4 <- _cons@a, ) ///
 (clsn5 <- _cons@a, ) ///
 (clsn6 <- _cons@a, ) ///
  , latent(Intercept ) ///
    cov(e.clsn1@b e.clsn3@b ///  // Constrain residual variances
        e.clsn4@b e.clsn5@b ///  // to be the same
        e.clsn6@b) ///
    nocapslatent
estimates store m1
di 5.27 / (5.27 + 4.68) // Random intercept model allows calculating an ICC
estat gof, stats(chi2 rmsea indices residuals)
                        // Poor model fit according to all tests
  // Likelihood ratio test:
lrtest m1 m0 // Massive improvement in fit for Model 1, though

// Model 2: Fixed intercept, fixed slope model (Table 2.4)
// Only mean intercept, intercept variance, and residual variance are 
// estimated



sem (Intercept@1 -> clsn1, ) ///  // Constrain paths to be 1
    (Intercept@1 -> clsn3, ) ///
    (Intercept@1 -> clsn4, ) ///
    (Intercept@1 -> clsn5, ) ///
    (Intercept@1 -> clsn6, ) ///
    (clsn1 <- _cons@a, ) ///     // Constrain intercepts to be the same
    (clsn3 <- _cons@a, ) ///
    (clsn4 <- _cons@a, ) ///
    (clsn5 <- _cons@a, ) ///
    (clsn6 <- _cons@a, ) ///
    (Slope@0 -> clsn1, ) ///     // Determine temporal structure
    (Slope@2 -> clsn3, ) ///
    (Slope@3 -> clsn4, ) ///
    (Slope@4 -> clsn5, ) ///
    (Slope@5 -> clsn6, ) ///
  , covstruct(_lexogenous, diagonal) ///
    latent(Intercept Slope ) ///
    cov(Intercept@0 ///        // Set intercept variance to 0
        Slope@0 ///            // Set slope variance to 0
 e.clsn1@b ///
        e.clsn3@b ///
        e.clsn4@b ///
        e.clsn5@b ///
        e.clsn6@b) ///
 means(Slope) ///           // Estimate slope
 nocapslatent
estimates store m2
estat gof, stats(chi2 rmsea indices residuals)
                     // Poor fit, even worse than Model 1

// Model 3: Random intercept, fixed slope (Table 2.5)



sem (Intercept@1 -> clsn1, ) ///  // Constrain paths to be 1
    (Intercept@1 -> clsn3, ) ///
    (Intercept@1 -> clsn4, ) ///
    (Intercept@1 -> clsn5, ) ///
    (Intercept@1 -> clsn6, ) ///
    (clsn1 <- _cons@a, ) ///     // Constrain intercepts to be the same
    (clsn3 <- _cons@a, ) ///
    (clsn4 <- _cons@a, ) ///
    (clsn5 <- _cons@a, ) ///
    (clsn6 <- _cons@a, ) ///
    (Slope@0 -> clsn1, ) ///     // Determine temporal structure
    (Slope@2 -> clsn3, ) ///
    (Slope@3 -> clsn4, ) ///
    (Slope@4 -> clsn5, ) ///
    (Slope@5 -> clsn6, ) ///
    , covstruct(_lexogenous, diagonal) ///
 latent(Intercept Slope ) ///
    cov(Slope@0 ///            // Set slope variance to 0
 e.clsn1@b ///
        e.clsn3@b ///
        e.clsn4@b ///
        e.clsn5@b ///
        e.clsn6@b) ///
 means(Slope) ///           // Estimate slope
 nocapslatent
estimates store m3
estat gof, stats(chi2 rmsea indices residuals)

  // Likelihood ratio test:
lrtest m3 m2 // Improvement in fit for Model 3

// Model 4: Random intercept, random slope (Table 2.6)


sem (Intercept@1 -> clsn1, ) ///  // Constrain paths to be 1
    (Intercept@1 -> clsn3, ) ///
    (Intercept@1 -> clsn4, ) ///
    (Intercept@1 -> clsn5, ) ///
    (Intercept@1 -> clsn6, ) ///
    (clsn1 <- _cons@a, ) ///     // Constrain intercepts to be the same
    (clsn3 <- _cons@a, ) ///
    (clsn4 <- _cons@a, ) ///
    (clsn5 <- _cons@a, ) ///
    (clsn6 <- _cons@a, ) ///
    (Slope@0 -> clsn1, ) ///     // Determine temporal structure
    (Slope@2 -> clsn3, ) ///
    (Slope@3 -> clsn4, ) ///
    (Slope@4 -> clsn5, ) ///
    (Slope@5 -> clsn6, ) ///
  , covstruct(_lexogenous, diagonal) ///
    latent(Intercept Slope ) ///
    cov(Intercept*Slope ///            // Include intercept-slope covariance
        e.clsn1@b ///
        e.clsn3@b ///
        e.clsn4@b ///
        e.clsn5@b ///
        e.clsn6@b) ///
 means(Slope) ///           // Estimate slope
 nocapslatent
estimates store m4
estat gof, stats(chi2 rmsea indices residuals)

  // Likelihood ratio test:
lrtest m4 m3 // Improvement in fit for Model 4

// Model 4 estimated as a multilevel model (Table 4.1):

// Read in data from clsn_cov.dat, this time using -corr2data-
#delimit ;
matrix input C = (
 6.3944,
 3.2716,  7.5282,
 4.1435,  6.0804, 10.7290,
 3.7058,  5.1597,  6.5672, 10.2920,
 4.1286,  5.7608,  7.2365,  7.6463, 12.9085,
-0.0940, -0.0390, -0.1521, -0.1104, -0.1469, 0.2502
 )
;
#delimit cr

corr2data clsn1 clsn3 clsn4 clsn5 clsn6 sex, n(851) ///
          means(37.9542 37.2785 37.0463 36.5696 36.1363 0.49) ///
    cov(C) cstorage(lower) clear
    
correlate, cov // Seems to work

gen id = _n // Create individual identifier

reshape long clsn, i(id) j(grade) // Convert to long format

mixed clsn grade || id: grade, var cov(uns)

Reference

Preacher, Kristopher J., Aaron L. Wichman, Robert C. MacCallum, and Nancy E. Briggs. 2008. Latent Growth Curve Modeling. Sage. doi: 10.4135/9781412984737