Multiple imputation of longitudinal data

This allows replicating the third example in Allison (2002, pp. 74-76).

// MI Example 3

use "https://statisticalhorizons.com/wp-content/uploads/hip.dta", clear
drop if wave == 4 // Not sure these are the correct data

xtset sid 

preserve
drop if missing(cesd, srh, adl, walk, pain)
bysort sid: drop if _N < 3
eststo clear
eststo: xtreg cesd srh walk adl pain ib3.wave, fe
restore

preserve

eststo: xtreg cesd srh walk adl pain ib3.wave, fe
restore

preserve
mi set mlong
mi register impute cesd srh walk adl pain wave

mi impute mvn cesd srh walk adl pain wave, ///
   add(10) burnin(500) burnbetween(30) 

eststo: mi estimate, post: xtreg cesd srh walk adl pain ib3.wave, fe
restore

preserve

reshape wide adl pain srh walk cesd, i(sid) j(wave)
mi set mlong
mi register impute cesd* srh* walk* adl* pain*

mi impute mvn cesd* srh* walk* adl* pain*, ///
   add(10) burnin(500) burnbetween(200)

mi reshape long adl pain srh walk cesd, i(sid) j(wave)
eststo: mi estimate, post: xtreg cesd srh walk adl pain ib3.wave, fe
restore

// Table 6.4

esttab, wide se nonumbers mtitle("LD by person" "LD by person-wave" "MI by person-wave" "MI by person")

Reference

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

Handling missing values in Stata

This allows replicating the analyses in Allison (2002, pp. 68-73).

// MI example 2
use spanking age educ income91 sex race marital region childs god using "C:\Users\User\Dropbox (FAMSIZEMATTERS)\methods and data\GSS1994.dta", clear
recode spanking (1 = 4) (2 = 3) (3 = 2) (4 = 1)
generate female  = (sex == 2)
generate black   = (race == 2)
recode income91 ( 1 =   500) ( 2 =  2000) ( 3 =  3500) ( 4 =  4500) ( 5 =  5500) /// 
                ( 6 =  6500) ( 7 =  7500) ( 8 =  9000) ( 9 = 11250) (10 = 13750) ///
                (11 = 16250) (12 = 18750) (13 = 21250) (14 = 23750) (15 = 27500) ///
                (16 = 32500) (17 = 37500) (18 = 45000) (19 = 55000) (20 = 67500) ///
                (21 = 75000), gen(income)
replace income = income / 1000

generate nochild = (childs == 0)         if !missing(childs)
generate nodoubt = (god == 6)            if !missing(god)
generate nevmar  = (marital == 5)        if !missing(marital)
generate divsep  = inlist(marital, 3, 4) if !missing(marital)
generate widow   = (marital == 2)        if !missing(marital)
generate east    = inlist(region, 1, 2)
generate midwest = inlist(region, 3, 4)
generate south   = inlist(region, 5, 6, 7)

misschk spanking female black income educ nodoubt nochild age east midwest south nevmar divsep widow

eststo clear
eststo: ologit spanking female black income educ nodoubt nochild age east midwest south nevmar divsep widow

drop if missing(marital)

preserve
recode educ (.d .n = .) 
recode spanking (.d .i .n = .)

mi set mlong
mi register imputed spanking female black income educ nodoubt nochild age east midwest south nevmar divsep widow

mi impute mvn spanking female black income educ nodoubt nochild age east midwest south nevmar divsep widow, ///
   add(5) burnin(500) burnbetween(200) emlog emoutput

foreach x of varlist female black nodoubt nochild east midwest south nevmar divsep widow {
   replace `x' = 0 if `x'  < .5 & _mi_m != 0
   replace `x' = 1 if `x' >= .5 & _mi_m != 0
}

replace spanking = 1 if                   spanking < 1.5 & _mi_m != 0
replace spanking = 2 if spanking >= 1.5 & spanking < 2.5 & _mi_m != 0
replace spanking = 3 if spanking >= 2.5 & spanking < 3.5 & _mi_m != 0
replace spanking = 4 if spanking >= 3.5                  & _mi_m != 0

eststo: mi estimate, post: ologit spanking female black income educ nodoubt nochild age east midwest south nevmar divsep widow
restore

preserve
recode educ (.d .n = .) 
recode spanking (.d .i .n = .)

mi set mlong
mi register imputed spanking female black income educ nodoubt nochild age east midwest south nevmar divsep widow

mi impute chained (mlogit) spanking (regress) income educ (logit) nodoubt nochild = female black age east midwest south nevmar divsep widow, ///
   add(5) burnin(20) force

eststo: mi estimate, post: ologit spanking female black income educ nodoubt nochild age east midwest south nevmar divsep widow
restore

preserve
mi set mlong
mi register imputed spanking female black income educ nodoubt nochild age east midwest south nevmar divsep widow

drop if missing(spanking)
mi impute chained (mlogit) spanking (regress) income educ (logit) nodoubt nochild = female black age east midwest south nevmar divsep widow, ///
   add(5) burnin(20) force

eststo: mi estimate, post: ologit spanking female black income educ nodoubt nochild age east midwest south nevmar divsep widow
restore

esttab, wide se keep(spanking:) nonumbers modelwidth(15) ///
  mtitle("Listwise deletion" "Normal data augmentation" "Sequential regression" "Seq. regression w/out missings") ///
  title(Coefficient estimates and standatd errors for cumulative logit models predicting SPANKING)

Reference

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

Using additional variables in multiple imputation

This allows replicating Table 6.2 in Allison (2002).

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

mi set mlong
mi register imputed csat act gradrat

eststo clear
eststo: regress csat


// Impute using ACT
mi impute mvn csat act, ///
   add(5) burnin(500) burnbetween(200) emlog emoutput

eststo: mi estimate, post: regress csat

// PCT25 is missing altogether in the data

// Impute using ACT and GRADRAT
mi impute mvn csat act gradrat, ///
   add(5) burnin(500) burnbetween(200) emlog emoutput

eststo: mi estimate, post: regress csat

esttab, not se mtitle("No imputation" "ACT" "ACT and GRADRAT") nonumbers ///
        coeflabel(_cons "Mean") modelwidth(15) title("Mean (and standard errors) of CSAT with different variables used in imputation")

Reference

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

Analysis of incomplete data using multiple imputation in Stata

This replicates MI Example 1 of Allison (2002, pp. 41-50) using Stata 14.

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

mi set mlong
mi register imputed gradrat csat private lenroll stufac rmbrd act
mi impute mvn gradrat csat private lenroll stufac rmbrd act, ///
   add(5) burnin(500) burnbetween(200) emlog emoutput ///
   saveptrace(trace, replace)

preserve
mi ptrace describe trace 
mi ptrace use trace, clear

// Generate regression coefficient
generate b = v_y2y1 / v_y2y2 

// Figure 5.1
twoway line b iter if inrange(iter, 1, 100), ///
       xtitle(Iteration) ytitle(b(csat)) ///
       ylabel(, format(%6.3f)) name(figure51, replace)


// Figure 5.2
tsset iter
ac v_y2y1, lags(100) ciopts(color(white)) note("") name(figure52, replace)

restore


// Table 5.3
eststo clear
foreach i of numlist 1/5 {
 qui eststo: regress gradrat csat lenroll private stufac rmbrd if _mi_m == `i'
}
esttab, not se wide nostar noobs order(_cons) varlabel(_cons "Intercept") nomtitle


// Figure 5.3
mi estimate: regress gradrat csat lenroll private stufac rmbrd

Reference

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

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

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

Dummy variable adjustment for missing values in Stata

This piece of code replicates Table 3.1 in Allison (2002).

clear
set seed 1

// Generate data
set obs 10000
drawnorm x z, ///
         corr(1, .5, 1) cstorage(lower) 
generate e = rnormal()
generate y = x + z + e

// Drop 1/2 of values from z
generate d = (runiform() > . 5)
generate zstar1 = z if d
replace  zstar1 = . if !d

// Substitute missing values
qui sum zstar1
generate zstar2 = zstar1
replace  zstar2 = r(mean) if !d

eststo clear
eststo: regress y x z
eststo: regress y x zstar1
eststo: regress y x zstar2 d

esttab using test.tex, b(2) not nostar nocons rename(zstar1 z zstar2 z) ///
        mtitles("Full data" "Listwise deletion" "Dummy variable adjustment") ///
        title(Regression in Simulated Data for Three Methods) replace booktabs 

Reference

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

Linear interpolation using ipolate

In some instances, it makes sense to fill in gaps in a dataset, e.g. a time-series, using linear interpolation. Stata offers the -ipolate- command for this. In a "long" data set containing the value value, where the variable country identifies the units i and the variable year the time scale j, the command:

by country, sort: ipolate value year, generate(interpolated_value)

creates a new variable interpolated_value that has gaps filled in with linearly interpolated values. However, any first or last observation of a time series will not be imputed, as this is considered to be extrapolation by the -ipolate- command. In order to have these time points filled in as well, the -epolate- option needs to be specified:

by country, sort: ipolate value year, generate(interpolated_value) epolate

A useful additional variable for sensitivity checks is a dummy indicator distinguishing between original and interpolated values:

// Create interpolation indicator
generate interpolated_value_dummy = (value == .)
label define interpolated_value_dummy 0 "Original value" 1 "Interpolated value"
label val interpolated_value_dummy interpolated_value_dummy
label var interpolated_value_dummy "Value interpolation y/n"

A general problem with these interpolated values is that when they are analyzed like complete data, standard errors will be underestimated (Allison 2002) even if the linear interpolations are consistent. An alternative approach (to which, however, the same caveat applies) would be to use Lowess smoothing (Cleveland 1979) via -lowess-.

References

Allison, Paul D. 2002. Missing Data. Sage.
Cleveland, William S. 1979. "Robust Locally Weighted Regression and Smoothing Scatterplots." Journal of the American Statistical Association 74(368):829-836. doi: 10.2307/2286407