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

Expectation Maximization (EM) for missing values using Stata

The code below allows replicating the analyses from Allison (2002, pp. 21-3).

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

// Table 4.1


eststo clear
estpost summarize gradrat csat lenroll private stufac rmbrd act

esttab using test.tex, cells("count(label(Nonmissing cases)) mean(label(Mean) fmt(2)) sd(label(SD) fmt(2))") ///
        nomtitle nonumber ///
        title(Descriptive Statistics for College Data Based on Available Cases) ///
        booktabs replace

// Tabe 4.2

eststo clear
eststo: regress gradrat csat lenroll private stufac rmbrd

#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 Listwise Deletion) 
         booktabs append ;
#delimit cr

// EM imputation
mi set mlong
mi register imputed gradrat csat lenroll private stufac rmbrd act 
mi impute mvn gradrat csat lenroll private stufac rmbrd act, emonly
matrix m = r(Beta_em)' // Transpose matrix of imputed means
matrix C = corr(r(Sigma_em)) // Matrix of correlations
matrix variances = diag((vecdiag(r(Sigma_em)))) // Matrix of variances
matrix sds = vecdiag(cholesky(variances))' // Vector of standard deviations
matrix descriptives = m, sds // Matrix needed for Table 4.3

// Table 4.3 

esttab matrix(descriptives, fmt(2 2)) using test.tex, ///
       nomtitle title("Means and Standard Deviations from the EM Algorithm") ///
       booktabs append

// Table 4.4

esttab matrix(C, fmt(3 3)) using test.tex, ///
       nomtitle title("Correlations from the EM Algorithm") ///
       booktabs append

// Table 4.5
drop *                                         // Get rid of data but not matrices
ssd init gradrat csat lenroll private stufac rmbrd act 
ssd set observations 1302
ssd set means (stata) m
ssd set sd (stata) sds
ssd set corr (stata) C

eststo clear

eststo: sem (gradrat <- csat lenroll private stufac rmbrd) 

#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 Based on the EM Algorithm)
  keep(gradrat:) eqlabels("", none) // Removes equation label
  booktabs append 
  ;
#delimit cr

Reference

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

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