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

Random graphs (69): ICC's with confidence intervals

use sharew1_rel2-6-0_gv_isced.dta, clear

// Reshape data into long format
reshape long iscedy_c, i(mergeid) j(child)
// Fix variable of interest
recode iscedy_c (-7 = .a "not yet coded (temporary)") ///
                (-2 = .b "refusal") ///
                (-1 = .c "don't know") ///
                (95 = .d "still in school") ///
                (97 = .e "other") ///
                ( . = .f "missing") ///
               , gen(years)
label var years "Years of education"
    
// Israel doesn't provide a ISCED-to-years conversion, thus it's dropped here
drop if country == 25

// Fix country variable
decode country, gen(cntry)

// Calculate ICC's per country
tempname foo
postfile `foo' str100 commandline_str str20 cntry icc icclb iccub N N_groups using "C:\Windows\Temp\test.dta", replace    
    
levelsof cntry, local(country)

foreach x of local country {
  *di "`x'"
  qui mixed years || mergeid: if cntry == "`x'", reml
  qui estat icc
  
  matrix groups = e(N_g)
  scalar n_g = groups[1,1]
  matrix cis =  r(ci2)
  scalar icclb = cis[1,1]
  scalar iccub = cis[1,2]
  
  post `foo' (e(cmdline)) ("`x'") (r(icc2)) (icclb) (iccub) (e(N)) (n_g)
  
  matrix stuff1 = (r(icc2), icclb, iccub, e(N), n_g)
  if "`x'" == "Austria" matrix table1 = stuff1
  else matrix table1 = (table1\stuff1)
}

postclose `foo'

matrix rownames table1 = `country'
matrix colnames table1 = "ICC" "CI lower" "CI upper" "Children" "Families"
esttab matrix(table1, fmt(2 2 2 0 0))

use "C:\Windows\Temp\test.dta", clear

egen order_ = rank(-icc), unique
labmask order_, value(cntry)

twoway (rcap   icc icc order_, horizontal) ///
       (rspike icclb iccub order_, horizontal) ///
      , ylabel(1/11, val) legend(off) ytitle("") xscale(alt) ///
        xtitle("Sibling correlations in educational attainment") ///
        note(" " "{it:Source:} SHARE wave 1, doi:10.6103/SHARE.w1.260", span)

erase "C:\Windows\Temp\test.dta"

Avoding the "initial values not feasible" (r(1400)) error with -melogit-

melogit home time || country: time, cov(uns) or

yields an error message "initial values not feasible." The same error occurs with the much simpler:

melogit home || country:, or

In order to further simplify estimation, the number of integration points can be reduced (the default is 7):

melogit home || country:, intpoints(2) or 

Drawing on the estimates of this model as starting values, the desired model can then be estimated:

melogit home || country:,  intpoints(2) or
mat a = e(b)
melogit home time || country:, intpoints(2) or from(a, skip)
mat b = e(b)
melogit home time || country: time, cov(uns) intpoints(2) or from(b, skip)

Random graphs (56): Intraclass Correlation Coefficients with confidence intervals

webuse productivity, clear

tempname foo
postfile `foo' str100 commandline_str str100 description_str icc icclb iccub N N_groups using "C:\Windows\Temp\test.dta", replace

qui mixed gsp || region: if year <= 1973
qui estat icc

matrix groups = e(N_g)
scalar n_g = groups[1,1]
matrix cis =  r(ci2)
scalar icclb = cis[1,1]
scalar iccub = cis[1,2]
local desc "ICC 1970 to 1973"

post `foo' (e(cmdline)) ("`desc'") (r(icc2)) (icclb) (iccub) (e(N)) (n_g)

qui mixed gsp || region: if year >= 1974 | year <= 1977
qui estat icc

matrix groups = e(N_g)
scalar n_g = groups[1,1]
matrix cis =  r(ci2)
scalar icclb = cis[1,1]
scalar iccub = cis[1,2]
local desc "ICC 1974 to 1977"

post `foo' (e(cmdline)) ("`desc'") (r(icc2)) (icclb) (iccub) (e(N)) (n_g)

qui mixed gsp || region: if year >= 1978 | year <= 1981
qui estat icc

matrix groups = e(N_g)
scalar n_g = groups[1,1]
matrix cis =  r(ci2)
scalar icclb = cis[1,1]
scalar iccub = cis[1,2]
local desc "ICC 1978 to 1981"

post `foo' (e(cmdline)) ("`desc'") (r(icc2)) (icclb) (iccub) (e(N)) (n_g)

qui mixed gsp || region: if year >= 1981 | year <= 1984
qui estat icc

matrix groups = e(N_g)
scalar n_g = groups[1,1]
matrix cis =  r(ci2)
scalar icclb = cis[1,1]
scalar iccub = cis[1,2]
local desc "ICC 1981 to 1984"

post `foo' (e(cmdline)) ("`desc'") (r(icc2)) (icclb) (iccub) (e(N)) (n_g)

qui mixed gsp || region: if year >= 1985 | year <= 1986
qui estat icc

matrix groups = e(N_g)
scalar n_g = groups[1,1]
matrix cis =  r(ci2)
scalar icclb = cis[1,1]
scalar iccub = cis[1,2]
local desc "ICC 1985 to 1986"

post `foo' (e(cmdline)) ("`desc'") (r(icc2)) (icclb) (iccub) (e(N)) (n_g)

postclose `foo'

use "C:\Windows\Temp\test.dta", clear
list

encode description_str, gen(description)
twoway (dot icc description, horizontal) ///
       (rcap icclb iccub description, horizontal)  ///
  , ylabel(1/5, val) ytitle("") legend(off) ///
    xlabel(0 (.1) 1) xtitle("ICC") name(figure, replace)
  
  
erase "C:\Windows\Temp\test.dta"

Firebaugh (1997): Analyzing Repeated Surveys in Stata (ch. 4)

Note how this is restricted to the linear decompositions!

Intracohort change versus overall change

use year race age cohort racmar using GSS7214_R4.DTA, clear

// Select sample
keep if inlist(year, 1974, 1984, 1994)
keep if race == 1   // Whites only
keep if age >= 20   
keep if !missing(age, racmar)

// Generate 10-year cohort variable
generate cohort10 = .
local lowerbracket = 1965
local upperbracket = 1974

foreach x of numlist 1/10 {
 replace cohort10 = `x' if cohort >= `lowerbracket' & cohort <= `upperbracket'
 label define cohort10 `x' "`lowerbracket' to `upperbracket'", modify
 local lowerbracket = `lowerbracket' - 10
 local upperbracket = `upperbracket' - 10
}

label value cohort10 cohort10
label var year "Year"

// Opposing interracial marriage
recode racmar (1 = 1 "Yes") (2 = 0 "No") (.d .i .n = .), gen(oppose)
label var oppose "Opposes interracial marriage"
generate oppose_perc = oppose * 100

// Create Table 4.1 without change
table cohort10 year, contents(mean oppose_perc n oppose_perc) format(%9.1f) center

// Different version of Table 4.1, including calculated changed
preserve
qui table cohort10 year, contents(mean oppose_perc n oppose_perc) format(%9.1f) center replace
rename table1 oppose_perc 
rename table2 oppose_n
reshape wide oppose_perc oppose_n, i(cohort10) j(year)
gen change1984 = oppose_perc1984 - oppose_perc1974 
gen change1994 = oppose_perc1994 - oppose_perc1984 
format change1984 change1994 %9.1f
list, noobs sep(0) abbreviate(20) subvar
restore

Empirical example for linear decomposition: Trend in antiblack prejudice

use year race cohort racdin racpush racseg racmar using GSS7214_R4.DTA, clear

// Create sample
keep if inlist(year, 1972, 1976, 1980, 1984)
keep if race == 1

// Create outcome
alpha racdin racpush racseg racmar, item
egen zracdin  = std(racdin)
egen zracpush = std(racpush)
egen zracseg  = std(racseg)
egen zracmar  = std(racmar)
scores y = total(zracdin zracpush zracseg zracmar), minvalid(3)
gen prejudice = 6 - y  // Difficult to figure out the exact outcome
drop y                 // But this seems to be close enough

// Calculate total change in prejudice
mean prejudice, over(year) coeflegend
local totalchange = _b[1984] - _b[1972]

// Calculate annual change within cohorts
reg prejudice year cohort
local b1 = _b[year]
local b2 = _b[cohort]

// Calculate contributions of intracohort change and cohort replacement
local intracohort = `b1' * (1984 -1972)
di `intracohort'

qui mean cohort, over(year) coeflegend
local cohortreplacement = `b2' * (_b[1984] - _b[1972])
di `cohortreplacement'

// Calculate estimated total change
local estim_change = `intracohort' + `cohortreplacement'
di `estim_change'

// Output results
matrix results = (`totalchange')\(`b1')\(`b2')\(`intracohort')\(`cohortreplacement')\(`estim_change')
esttab matrix(results, fmt(%9.2f)), coeflabel(r1 "Total change" ///
                                              r2 "Intracohort slope" ///
                                              r3 "Intercohort slope" ///
                                              r4 "Estimated contribution of intracohort change" ///
                                              r5 "Estimated contribution of cohort replacement" ///
                                              r6 "Estimated total change") ///
                                              mtitle("Result")  ///
                                              varwidth(45)

Empirical example of the same-sign rule: Gender role attitudes

use year cohort fework fepres fepol fehome if year >= 1972 & year <= 1988 using GSS7214_R4.DTA, clear

// Prepare variables
recode fework fepres (2 = 0)
recode fepol fehome  (2 = 1) (1 = 0)

// Create sample for Table 4.2
preserve
keep if inlist(year, 1972, 1974, 1988)
fre year
gen id = _n
reshape wide fework fepres fepol fehome, i(id) j(year)
replace fepol1972 = fepol1974
replace fehome1972 = fehome1974
reshape long
drop if year == 1974
recode year (1972 = 0 "1972") (1988 = 1 "1988"), gen(yr)

// Table 4.2

foreach y of varlist fework fepres fepol fehome {
 qui regress `y' yr
 local mean1972 = _b[_cons]
 local mean1988 = _b[_cons] + _b[yr]
 local differen = _b[yr]
 local t        = _b[yr] / _se[yr]
 *di "`y'" _skip(5) `mean1972' _skip(5) `mean1988' _skip(5) `differen' _skip(5) `t'
 matrix stuff1 = (`mean1972' , `mean1988' , `differen' , `t')
 if "`y'" == "fework" matrix table42 = stuff1
 else matrix table42 = (table42\stuff1)
}
 
matrix rownames table42 = "WORK" "PRES" "POLI" "HOME"
matrix colnames table42 = "Mean 1972" "Mean 1988" "Change" "t-value"
esttab matrix(table42, fmt(3 3 3 1))

restore

// Table 4.3


foreach y of varlist fework fepres fepol fehome {
 qui logit `y' year cohort
 local n             = e(N)
 local intracohort   = _b[year]
 local intracohort_t = _b[year] / _se[year]
 local intercohort   = _b[cohort]
 local intercohort_t = _b[cohort] / _se[cohort]
 *di "`y'" _skip(5) `n' _skip(5) `mean1988' _skip(5) `differen' _skip(5) `t'
 matrix stuff = (`n', `intracohort', `intracohort_t', ///
                      `intercohort', `intercohort_t')
 if "`y'" == "fework" matrix table43 = stuff
 else matrix table43 = (table43\stuff)
}

matrix rownames table43 = "WORK" "PRES" "POLI" "HOME"
matrix colnames table43 = "N" "Within cohort" "t" "Cross-cohort" "t"

esttab matrix(table43, fmt(0 %9.3f %9.1f %9.3f %9.1f)), modelwidth(14)


Reference

Firebaugh, Glenn. 1997. Analyzing Repeated Surveys. Sage. doi: 10.4135/9781412983396