Instrumental variable regression

clear
set seed 1
set obs 10000
eststo clear

// Simulate data for IV regression
matrix input c = ( 1, .5, 0, .5\ ///
                  .5,  1, 0,  0\ ///
                   0,  0, 1,  0\ ///
                  .5,  0, 0,  1)
corr2data x u e z, corr(c)

// Calculate y
generate y = 1 + 2*x + 2*u + 5*e

// Without omitted variable bias
eststo: regress y x u

// With omitted variable bias
eststo: regress y x

// Reduced form
eststo: regress y z
local rf = _b[z]

// First stage
eststo: regress x z
local fs = _b[z]
predict xhat
predict xres, resid

// Two-stage predictor substitution (2SPS)
eststo: regress y xhat

// Two-stage residual inclusion (2SRI)
eststo: regress y x xres

// Two-stage least square (2SLS)
eststo: ivregress 2sls y (x = z)

// IV estimate = reduced form/first stage
di `rf'/`fs'



// Table
esttab, b(2) se(2) mtitle("DGP" "Bias" "Reduced form" "1st stage" "2SPS" "2SRI" "2SLS") ///
        coeflabel(x "Predictor x" ///
                  z "Instrument z" ///
                  xres "1st stage residual" ///
                  _cons "Intercept") ///
        rename(xhat x) drop(u) ///
        stats(N, label(Observations) fmt(%9.0gc)) ///
        title(Different varieties of instrumental variable regression) ///
        varwidth(20)

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

Replicating simple analyses with -corr2data-

A correlation table is a pivotal feature of academic papers in many disciplines such as psychology. The reason for this is that many basic statistical methods rely on correlation matrices only. Thus, the vector of means and standard deviations plus the correlation matrix and the sample size are sufficient to replicate many standard analyses such as OLS regression or principal component analysis. -corr2data- allows to do so by creating an artficial data set based on the correlation matrix (or, covariance matrix), means, standard deviations, and sample size, which can then be used for replicating or modifying the analyses.

Table 1 of Huff-Corzine et al. (1986)

In order to try out -corr2data-, we will replicate a randomly chosen research paper, namely Huff-Corzine et al. (1986). Huff-Corzine and colleagues deal with the North–South differences in homicide in the US, and, more importantly for our purposes, they provide a full correlation table and conduct OLS regression analyses.


First, we need to read in the correlation matrix as reported in Figure 1. (N.B. Stata 12 now comes with the -ssd- command which might allow for a more convenient way of reading in summary statistics.)

version 12
clear
capture which estout
if _rc ssc install estout

// Read in correlation table
#delimit ;
input str30 varnames double (a1-a9);
"1 Homicide rate"
    1.00  0.91  0.83  0.88  0.43 -0.06 -0.30  0.65  0.83;
"2 Structural poverty index"
    0.91  1.00  0.78  0.90  0.40 -0.04 -0.17  0.65  0.84;
"3 Southernness index"
    0.83  0.78  1.00  0.72  0.34  0.14 -0.44  0.71  0.93;
"4 Perc nonwhite"
    0.88  0.90  0.72  1.00  0.32  0.08 -0.14  0.68  0.76;
"5 Perc ages 20-34"
    0.43  0.40  0.34  0.32  1.00 -0.14 -0.21 -0.09  0.32;
"6 Perc rural"
   -0.06 -0.04  0.14  0.08 -0.14  1.00  0.14  0.28  0.10;
"7 Hospital beds/100K"
   -0.30 -0.17 -0.44 -0.14 -0.21  0.14  1.00 -0.20 -0.30;
"8 Gini index"
    0.65  0.65  0.71  0.68 -0.09  0.28 -0.20  1.00  0.73;
"9 Perc born in South"
    0.83  0.84  0.93  0.76  0.32  0.10 -0.30  0.73  1.00;
end;
#delimit cr

// Convert to a matrix called M
mkmat a1-a9, matrix(M) 
matrix list M

// Read in variable names
levelsof(varnames), local(names)

local ednames = ""
  // Create local "ednames" that includes all variable names,
  // use function -strtoname()- to bring variable names into 
  // proper Stata format
foreach x of local names {
  *di strtoname("`x'", 1)
  local ed_name1 = strtoname("`x'", 1)
  *di "`ed_name1 '"
  local ed_name2 = substr("`ed_name1'", 4, .)
  *di "`ed_name2 '"
  local ednames = "`ednames'" + " " + "`ed_name2'"
}

In this step, we create the artificial data set based on the correlation matrix and means and s.d.'s:

// Create data based on correlation matrix M, vectors of
// means and sd's, specify no. of cases and variable names
corr2data `ednames', n(48) clear ///
     corr(M) ///
    means(7.12 7.12 17.71 10.42 19.75 34.22 7.81 .38 7.18) ///
   sds(4.25 3.88  9.22  8.85  1.27 14.37 1.57 .02 1.44)

One shouldn't underestimate the importance of fine-looking labels in a data set:

// Create nice-looking value labels based on the variable names
local varlabels = ""

foreach var of varlist _all {
  local varname = "`var'"
  *di "`varname'"
  local varname = subinstr("`varname'", "_"   , " ", .)
  *di "`varname'"
  local varname = subinstr("`varname'", "Perc", "%", .)
  *di "`varname'"
  lab variable `var'  "`varname'"
  local varlabels = "`varlabels'" + " " + "`varname'"
}

A brief check shows that every thing seems to be alright:

corr _all

tabstat _all, statistics(mean sd)

However, replicating the analyses yields different results:

eststo clear

// Run regression

eststo: regress Homicide_rate ///
          Structural_poverty_index Southernness_index ///
    Perc_nonwhite Perc_ages_20_34 Perc_rural ///
    Hospital_beds_100K Gini_index, beta  
    
eststo: regress Homicide_rate ///
          Structural_poverty_index ///
    Perc_nonwhite Perc_ages_20_34 Perc_rural ///
    Hospital_beds_100K Gini_index Perc_born_in_South, beta
    
estadd beta: est1 est2

esttab, r2 nonumbers mlabels("South index" "% born South") ///
        cells ((b(fmt(3) star) beta(fmt(3)))) label not)

Download the do-file here.

These results differ substantially from those reported by Huff-Corzine and others (see below), however I have a feeling that I should rather be trusting my own result.

Table 2 of Huff-Corzine et al. (1986)

A robustness check offered by -corr2data- is the -seed(#)- option, which allows for generating different artificial data sets. Specifying different seeds and comparing the results can serve as a check whether the summary statistics are really sufficient for replicating an analysis. Using different seeds in our example always leads to the same result.

Reference

Huff-Corzine, Lin, Jay Corzine, and David C. Moore. 1986. "Southern Exposure. Deciphering the South's Influence on Homicide Rates." Social Forces 64(4):906-924. doi: 10.1093/sf/64.4.906