Decomposing the difference between two means

// GSS 1990-2004
use year race sibs reg16 educ maeduc if inrange(year, 1990, 2004) ///
  & year != 2002 ///  // 2002 measures race in a non-standard way
    using "gss7221_r1.dta", clear

// In Table 7.8 Treiman states that he's got 17,090 cases (14,985 non-black + 
// 2,105 black). There seems to be no way I can get there with the data
// because maeduc, mother's education, only has 15,996 valid observations.

// The reason for this seems to be that in his do-file, he starts with
// the 1990 GSS file and then appends all other GSS year files -- including
// the 1990 file, thus including the 1990 file twice. Running the do-file with

// expand 2 if year == 1990

// uncommented (almost) replicates the numbers reported in the book.

drop year // Not needed

// Race variables
generate black    = (race == 2)
generate nonblack = !black
drop race

// Truncate number of siblings at 15
replace sibs = 15 if sibs > 15  & !missing(sibs)
label var sibs "Sibsize"

// South
gen south = (inrange(reg16, 5, 7))
label var south "Southern origin"
drop reg16

// Education
label var educ "Education"
label var maeduc "Mother's education"

// Listwise deletion
mark touse if !missing(educ, maeduc, sibs)
keep if touse



// Table 7.8a
eststo clear
local vlist educ maeduc sibs south 
local upper

local lower `vlist'

foreach v of local vlist {
  estpost correlate `v' `lower' if !black
  local nnonblack = e(N)
  foreach m in b rho p count {
   matrix `m' = e(`m')
  }
  if "`upper'"!="" {
   estpost correlate `v' `upper' if black
   local nblack = e(N)
   foreach m in b rho p count {
    matrix `m' = e(`m'), `m'
   }
  }
  ereturn post b
  foreach m in rho p count {
   quietly estadd matrix `m' = `m'
  }
  eststo `v'
  local lower: list lower - v
  local upper `upper' `v'
}

// Table 7.8b
esttab using table2.tex, replace nonumbers noobs mtitles not booktabs label nostar ///
       title(Correlations between study variables: Blacks (\emph{N} = `nblack') above, non-blacks (\emph{N}  = `nnonblack') below the diagonal) ///
mtitle("Education" "Mother's education" "Sibsize" "Southern origin")

eststo clear
eststo: estpost sum educ maeduc sibs south if black
eststo: estpost sum educ maeduc sibs south if !black

esttab using table1.tex, booktabs replace ///
       cells("mean(label(Mean) fmt(2))" "sd(label(SD) fmt(2) par)") ///
       label mtitle("Blacks" "Non-Blacks") ///
       title("Means and standard deviations of study variables")



// Table 7.9 
eststo clear
eststo: regress educ maeduc sibs south if touse & black
eststo: regress educ maeduc sibs south if touse & !black

esttab using table3.tex, booktabs replace ///
       b(2) se(2) r2(2) ///
       label mtitle("Blacks" "Non-Blacks") ///
       title("Coefficients of a model of years of schooling, for blacks and non-blacks, US adults, 1990--2004")



// Table 7.10
eststo clear
eststo: oaxaca educ maeduc sibs south if touse, by(black) detail noisily
eststo: oaxaca educ maeduc sibs south if touse, by(nonblack) detail noisily
esttab using table4.tex, booktabs replace ///
       b(2) nose not label mtitle("Blacks as reference" "Non-Blacks as reference") ///
       eqlabel("\emph{Overall}" ///
               "\emph{Differences in assets}" ///
               "\emph{Differences in returns to assets}" ///
               "\emph{Interactions}") ///
       coeflabels(overall:difference "Difference in years of schooling" ///
                 overall:endowments "Total due to difference in assets" ///
                 overall:coefficients "Total due to difference in returns" ///
                 overall:interaction "Total due to interactions") ///
       drop(overall:group*) ///
       varwidth(30) alignment(D{.}{.}{-1}) ///
       title("Decomposition of the difference in the mean years of schooling by blacks and non-blacks, US adults, 1990--2004")

Inverse probability of treatment weighting

// Allbus 2018
use xr19 sex age isced97 eastwest german pt09 pt10 wghtpew bik using "ZA5272_v1-0-0.dta", clear

// Outcome: Trust in media
recode pt09 pt10 (-9 = .)
alpha pt09 pt10 
scores trust = mean(pt09 pt10), nv(1)
label var trust "Trust in media"
drop pt09 pt10 

// Treatment: Internet use
gen internet = (xr19 == 1) if !inlist(xr19, -9, -8)
label define internet 1 "Internet user" 0 "Internet non-user"
label val internet internet
drop xr19

// Gender
gen female = (sex == 2)
label var female "Female sex"
drop sex

// Age
recode age (-32 = .)
label var age "Age"

// Education
recode isced97 (-32 = .)
label define isced97 1 "Basic" 2 "Lower secondary" 3 "Upper secondary" 4 "Post-secondary" 5 "Tertiary first" 6 "Tertiary second", modify
label var isced97 "Education"

// East
generate east = (eastwest == 2)
label var east "Eastern Germany"
drop eastwest

// German
recode german (1 = 0) (2 = 0) (3 = 1) (-50 = 1), gen(foreign)
label var foreign "Foreign citizen"
drop german

// Region
recode bik (-34 = .)
label define bik  1 "-1,999 inh." ///
                  2 "2,000-4,999 inh." ///
                  3 "5,000-19,999 inh." ///
                  4 "Zone 1-4, -50,000 inh." ///
                  5 "Zone 2-4, -100,000 inh." ///
                  6 "Zone 1, -100,000 inh." ///
                  7 "Zone 2-4, -500,000 inh." ///
                  8 "Zone 1, -500,000 inh." ///
                  9 "Zone 2-4, 499,999+ inh." ///
                 10 "Zone 1, 499,999+ inh.", modify
label var bik "Commuting zone"

// Keep only complete cases
keep if !missing(trust, internet, female, age, isced97, east, foreign, bik)

// See how imbalanced covariates are
  // Generate dummies for first table
foreach x of varlist isced97 bik {               // Plug in categorical variables here
  qui tab `x', gen(`x'x)                         // Create dummy variables
  foreach var of varlist `x'x* {
   local lab `: var label `var''
   *di "`lab'"                                   // Display label
   *di strpos("`lab'", "==") 
   local i = strpos("`lab'", "==") + 2
   local lab `: di substr("`lab'",  `i', .)'
   di "`lab'"
   label var `var' "  `lab'"                // Add space to label
  }
}
  // Calculate means
eststo clear
estpost tabstat trust female isced97x* age foreign east bikx*, by(internet) ///
                                                       statistics(mean sd) ///
                                                       columns(statistics) ///
                                                       casewise
   // Look at table
esttab, unstack cells(mean(fmt(2)) sd(fmt(2) par keep(trust age))) nonumber ///
        label collabels("Mean (SD)/Prop.") modelwidth(25) varwidth(25) ///
        refcat(isced97x1 "Education:" bikx1 "Communting zone:", nol) 
drop isced97x* bikx* // Drop variables created for table
eststo clear

// Alternative approach
  // Run desired model first
qui teffects ipw (trust) (internet east##(c.age##c.age) female foreign ib3.isced97 i.bik, logit)
  // Table of covariates
tebalance summarize, baseline // Same table as created above
matrix before = r(table)
matrix list before
  // Table of covariates after weighting
tebalance summarize
matrix after = r(table)
matrix list after
  // Merge tables
matrix both = before, after
matrix list both

  // Look at table
esttab matrix(both, fmt(2)), ///
       label nomtitle ///
       coeflabel(1.east "Eastern Germany" ///
                 age "Age" ///
                 c.age#c.age "Age X age" ///
                 1.east#c.age "Eastern Germany X age" ///
                 1.east#c.age#c.age "Eastern Germany X age X age") ///
       refcat(1.isced97 "Education:" 2.bik "Commuting zone:", nol) ///
       varwidth(30) modelwidth(18) ///
       collabel("Control (mean)" "Treated (mean)" "Control (var.)" "Treated (var.)" "Std. diff." "Std. diff. weighted" "Var. ratio" "Var. ratio weighted")

// Omnibus test: tests hypothesis that treatment model balances the covariates
tebalance overid

// Check overlap
teffects overlap, xtitle(Propensity score Internet non-user) ///
                  ytitle(Density) ///
                  legend(order(1 "Internet non-user" 2 "Internet user") pos(2) ring(0)) 

// IPTW
eststo: teffects ipw (trust) (internet east##(c.age##c.age) female foreign ib3.isced97 i.bik, logit), ate
eststo: teffects ipw (trust) (internet east##(c.age##c.age) female foreign ib3.isced97 i.bik, logit), atet

// Doing things by hand:
// Calculate weight for IPTW, ATE
logit internet east##(c.age##c.age) female foreign ib3.isced97 i.bik
predict probab
generate wate = .
replace  wate = 1/probab       if internet == 1
replace  wate = 1/(1 - probab) if internet == 0

// Calculate weight for IPTW, ATT
generate watt = .
replace  watt = 1                   if internet == 1
replace  watt = probab/(1 - probab) if internet == 0

// IPTW by hand
eststo: regress trust i.internet [pw = wate]
eststo: regress trust i.internet [pw = watt]
esttab, b(2) se(2) keep(main:) drop(_cons) rename(r1vs0.internet 1.internet) ///
        nobaselevels mtitles("IPTW ATE" "IPTW ATT" "IPTW ATE" "IPTW ATT") ///
        mgroup("teffects" "regress", pattern(1 0 1 0))

// Doubly robust
eststo: teffects ipwra (trust          east##(c.age##c.age) female foreign ib3.isced97 i.bik) ///
                     (internet east##(c.age##c.age) female foreign ib3.isced97 i.bik, logit)
eststo: teffects ipwra (trust          east##(c.age##c.age) female foreign ib3.isced97 i.bik) ///
                     (internet east##(c.age##c.age) female foreign ib3.isced97 i.bik, logit), atet
esttab, b(2) se(2) keep(main:) drop(_cons) rename(r1vs0.internet 1.internet) ///
        nobaselevels mtitles("IPTW ATE" "IPTW ATT" "IPTW ATE" "IPTW ATT" "IPTW RA ATE" "IPTW RA ATT") ///
        mgroup("teffects ipw" "regress" "teffects ipwra", pattern(1 0 1 0 1 0))

// Doubly robust by hand
tempvar internet0a internet1a internet0b internet1b ate att

// ATE
regress trust east##(c.age##c.age) female foreign ib3.isced97 i.bik if !internet [pw = wate]
predict `internet0a'
regress trust east##(c.age##c.age) female foreign ib3.isced97 i.bik if internet [pw = wate]
predict `internet1a'
generate `ate' = (`internet1a' - `internet0a')
summarize `ate'

// ATT
regress trust east##(c.age##c.age) female foreign ib3.isced97 i.bik if !internet [pw = watt]
predict `internet0b'
regress trust east##(c.age##c.age) female foreign ib3.isced97 i.bik if internet [pw = watt]
predict `internet1b'
generate `att' = (`internet1b' - `internet0b')
summarize `att'

Descriptive statistics table using -esttab-

// Open Allbus 1984-2016
use age iscd975 isei* sex using "ZA4586_v1-0-0.dta", clear

// Age
recode age (-32 = .)
label var age "Age"

// Education
recode iscd975 (-32 = .) ///
               (1 2 = 0 "Low") ///
               (3 4 = 1 "Medium") ///
               (5 = 2 "High"), gen(educ)
label var educ "Education"

// ISEI
generate isei = isei88  if inrange(isei88, 16, 90)
replace  isei = isei08  if missing(isei) & inrange(isei08, 16, 90)
replace  isei = isei68  if missing(isei) & inrange(isei68, 16, 90)
replace  isei = isei88a if missing(isei) & inrange(isei88a, 16, 90)
replace  isei = isei08a if missing(isei) & inrange(isei08a, 16, 90)
replace  isei = isei68a if missing(isei) & inrange(isei68a, 16, 90)
label var isei "ISEI" 

// Sex
recode sex (1 = 0 "Male") (2 = 1 "Female"), gen(female)
label var female "Sex"

// Listwise deletion
capture drop touse
mark touse
markout touse age educ sex isei

// Create dummies of categorical variables for descriptives table
foreach x of varlist educ female {               // Plug in categorical variables here
  qui tab `x', gen(`x'x)                         // Create dummy variables
  foreach var of varlist `x'x* {
   local lab `: var label `var''
   *di "`lab'"                                   // Display label
   *di strpos("`lab'", "==") 
   local i = strpos("`lab'", "==") + 2
   local lab `: di substr("`lab'",  `i', .)'
   label var `var' "\dumm `lab'"                // Add \dumm to label
  }
}

// Add this to Latex code to indent dummy variables
// % Indent for tables
// \newcommand*{\dumm}{\hspace*{0.5cm}}%

// Calculate descriptives and save them internally
eststo clear
estpost tabstat isei educx* age femalex* if touse, by(female) ///
                                                   statistics(mean sd min max) ///
                                                   columns(statistics)

// Format number of cases for table caption
local n: display %9.0gc e(N)
// Create table
esttab using "table1.tex", cells("mean(fmt(2) label(Prop./Mean)) sd(fmt(2) label(\emph{SD}) keep(isei age))") ///
                         refcat(educx1 "Education:" ///
                                femalex1 "Sex:", nolabel) ///
                         coeflabel(isei "Occupational status (ISEI)") ///
                         unstack ///
                         order(Total:* Male:* ) ///
                         nogap eqlabels(, span prefix(\multicolumn{@span}{c}{) suffix(}) erepeat(\cline{@span})) nonumber replace label ///
                         title("Descriptive statistics, \emph{N}~=~`n' \label{tab1}") ///
                         stats(N, fmt(%9.0gc) label(\emph{N})) ///
                         booktabs 
drop educx* femalex*  // Remove dummy variables created for table

Chapter 2 of Singer and Willett's (2003) book on longitudinal data analysis

There is a different take on the chapter here, but I like mine better.

// Open data set
use http://www.ats.ucla.edu/stat/stata/examples/alda/data/tolerance, clear

// Figure 2.1
list, sep(0) noobs
reshape long tol, i(id) j(age)
list if inlist(id, 9, 45, 1653), sepby(id) noobs

// Table 2.1
reshape wide
quietly estpost cor tol*, matrix listwise
esttab, unstack not noobs compress nostar b(%6.2f)

// Figure 2.2

reshape long
twoway (scatter tol age), by(id) ytitle(Tolerance) xtitle(Age) name(Figure22, replace)

// Figure 2.3

twoway (scatter tol age) ///
       (lowess tol age) ///
      , by(id, note("Graphs by id, lowess curves") legend(off)) ///
        ytitle(Tolerance) xtitle(Age) name(Figure23, replace)

// Table 2.2
preserve
tempfile table22
generate time = age - 11
statsby initial    = _b[_cons]  ///
        initial_se = _se[_cons] ///
        change       = _b[time] ///
        change_se    = _se[time] ///
        residual   = (e(rmse)^2) ///
        explained  = e(r2) ///
      , by(id) saving(`table22', replace): regress tol time
drop time
reshape wide  

merge 1:1 id using `table22', nogenerate
format    initial initial_se change change_se residual explained      exposure %6.2f
keep   id initial initial_se change change_se residual explained male exposure
order  id initial initial_se change change_se residual explained male exposure
list, sep(0) noobs abbrev(15)

tempfile figure28
save `figure28', replace

// Figure 2.4
stem initial, round(.01)
stem change, round(.01)
stem residual, round(.01)
stem explained, round(.01)
restore

// Figure 2.5

twoway (scatter tol age) ///
       (lfit tol age) ///
      , by(id, note("Graphs by id, OLS curves") legend(off)) ///
        ytitle(Tolerance) xtitle(Age) name(Figure25, replace)
   
// Figure 2.6

xtset id age
xtline tol, overlay t(age) i(id) legend(off) ///
            ytitle(Tolerance) ylabel(1 (1) 4) ///
            addplot(lowess tol age, lwidth(thick) lpattern(solid)) ///
            xtitle(Age) title("Observed data and lowess smoother") ///
            xsize(3) name(Figure26A, replace) nodraw
quietly regress tol i.id##c.age
predict tolhat
xtline tolhat, overlay t(age) i(id) legend(off) ///
               ytitle(Tolerance) ylabel(1 (1) 4) ///
               addplot(lfit tol age, lwidth(thick) lpattern(solid)) ///
               xtitle(Age) title("OLS trajectories") ///
               xsize(3) name(Figure26B, replace) nodraw
graph combine Figure26A Figure26B, col(2) name(Figure26, replace)

// Table 2.2
preserve
use `table22', clear
tabstat initial change, stat(mean sd) format(%6.2f)
cor initial change
restore

// Figure 2.7

quietly sum exposure, detail
generate highexposure = (exposure >= r(p50))
generate  lowexposure = (exposure  < r(p50))
generate  female      = !male
label var male         "Males"
label var female       "Females"
label var lowexposure  "Low exposure"
label var highexposure "High exposure" 

foreach z of varlist female male lowexposure highexposure {
  capture drop tolhat
  quietly regress tol i.id##c.age if `z'
  predict tolhat if `z'
  xtline tolhat if `z', ///
         overlay t(age) i(id) legend(off) ///
         ytitle("Predicted tolerance") ylabel(1 (1) 4) ///
         addplot(lfit tol age if `z', lwidth(thick) lpattern(solid)) ///
         xtitle(Age) title(`: variable label `z'') ///
         xsize(3) ysize(2) name(Figure27_`z', replace) nodraw
}
graph combine Figure27_female Figure27_male Figure27_lowexposure Figure27_highexposure, col(2) name(Figure27, replace)

// Figure 2.8

use `figure28', clear

quietly correlate initial male 
local pmca = round(r(rho), .01)
quietly correlate initial exposure 
local pmcb = round(r(rho), .01)
quietly correlate change male 
local pmcc = round(r(rho), .01)
quietly correlate change exposure 
local pmcd = round(r(rho), .01)

label define male -1 " " 1 "Male" 0 "Female" 2 " ", modify
label val male male

twoway dot initial male, ///
       ytitle("Predicted intercept") ///
       xtitle("Gender") ///
       xlabel(-1 0 1 2, val) ///
       ysize(2) xsize(3) ///
       note("{it:r} = `pmca'") ///
       name(Figure28A, replace) nodraw
twoway scatter initial exposure, ///
       ytitle("Predicted intercept") ///
       xtitle("Exposure") ///
       ysize(2) xsize(3) ///
       note("{it:r} = `pmcb'") ///
       name(Figure28B, replace) nodraw
twoway dot  change  male, ///
       ytitle("Predicted change") ///
       xtitle("Gender") ///
       ysize(2) xsize(3) ///
       xlabel(-1 0 1 2, val) ///
       note("{it:r} = `pmcc'") ///
       name(Figure28C, replace) nodraw
twoway scatter change  exposure, ///
       ytitle("Predicted change") ///
       xtitle("Exposure") ///       
       ysize(2) xsize(3) ///
       note("{it:r} = `pmcd'") ///    
       name(Figure28D, replace) nodraw
graph combine Figure28A Figure28B Figure28C Figure28D, col(2) row(2) name(Figure28, replace)

Reference

Singer, Judith D., and John B. Willett. 2003. Applied Longitudinal Data Analysis. Modeling Change and Event Occurrence. Oxford University Press. doi: 10.1093/acprof:oso/9780195152968.001.0001

Latex tables using -esttab-

// Table 1
eststo clear
bysort male country: eststo: estpost summarize hwhrs phwhrs gap ///
                                               workhours pworkhours ///
                                               genderroles femaleemployment ///
                                               married higherincome phigherincome ///
                                               age preschool schoolage ///
                                               retired pretired, listwise

esttab using test.tex, main(mean) aux(sd) label nodepvar nostar nonote nonumbers ///
                       b(%9.2f) varwidth(35) compress booktabs ///
                       title("Descriptive statistics: Means and standard deviations by gender and country") ///
                       mgroups("Women" "Men", pattern(1 0 0 1 0 0) ///
                       prefix(\multicolumn{@span}{c}{) suffix(})   ///
                       span erepeat(\cmidrule(lr){@span})) replace 

// Table 2
eststo clear
levelsof country, local(country)
foreach i of local country {
  capture drop touse
  mark touse
  markout touse   hwhrs workhours pworkhours phwhrs ///
                  genderroles femaleemployment ///
                  i.relativeincome pretired retired ///
                  age preschool schoolage married

  eststo: regress hwhrs workhours pworkhours phwhrs ///
          if female == 1 & country == `i' & touse
  eststo: regress hwhrs workhours pworkhours phwhrs ///
                  genderroles femaleemployment ///
                  i.relativeincome pretired retired ///
                  age preschool schoolage married ///
                  if female == 1 & country == `i' & touse
}

esttab using test.tex, se ar2 label varwidth(35) b(%9.2f) compress nodepvars order(_cons) booktabs ///
        mtitles("Germany" "Germany" "USA" "USA" "Finland" "Finland") ///
        title("Determinants of weekly housework hours of women in West Germany, USA, and Finland. OLS coefficients") ///
        append

Publication-style correlation tables in Stata

Quite a number of Stata users have engaged in programming commands to create readymade correlation tables to use in publications. Stata's own -correlate- and -pwcorr- lack many desirable features; for instance, means and standard deviations cannot be included in the table automatically and columns cannot be numbered. Furthermore, variable labels cannot be displayed, only variable names can be shown in the table.

sysuse auto
correlate price-foreign, means

-estpost-, -mkcorr-, -corrtab-, and -makematrix- are all user-written commands which aim to improve the Stata default table.

 

-estpost-

-estpost- stems from Ben Jann's mighty -estout- package. The general use of the command is as follows:

sysuse auto
capture which estout
if _rc ssc install estout
estpost correlate price-trunk, matrix
esttab, unstack not nonum compress noobs

In order to display both the upper and the lower triangle of the matrix, the -nohalf- option of -estpost- needs to be used:

estpost correlate price-trunk, matrix nohalf

However, getting -esttab- to use variable labels or to include descriptive statistics doesn't seem to be a trivial task, although the latter can be done. Finally, it should be noted that the default of -estpost correlate- is pairwise deletion, not listwise deletion; so that it is -pwcorr- rather than -correlate-. In order to obtain a correlation matrix based on listwise deletion, the -listwise- option needs to be specified.

 

-mkcorr-

-mkcorr- by Glenn Hoetker comes with many nice features and only few drawbacks.

sysuse auto, clear
capture which mkcorr
if _rc ssc install mkcorr
mkcorr price-trunk, log(corr_table)

It allows including summary statistics (mean, SD, minimum, and maximum) into the table:

mkcorr price-trunk, log(corr_table) replace means

The use of variable labels:

mkcorr price-trunk, log(corr_table) replace lab

The use of numbers in the column headers:

mkcorr price-trunk, log(corr_table) replace num

The inclusion of p-values:

mkcorr price-trunk, log(corr_table) replace means sig 

Furthermore, the number of decimal places for correlation can be manipulated via -cdec(#)- (-mdec(#)- for summary statistics).
The main drawback is that the table is not being displayed in the results window; instead it can only be written into a tab-separated text-only file that then needs additional formatting in Word or Excel. Also, the p-values are rather ugly; an option to include stars to denote statistical significance would be a great improvement. Again, the default setting is to provide correlations based on pairwise deletion; the option -casewise- will yield results based on listwise deletion.

-corrtab-

-corrtab- by Fred Wolfe is a tool that has a somewhat more limited functionality, or, more functions that I don't find to be of great use.

sysuse auto, clear
capture which corrtab
if _rc ssc install corrtab
corrtab price-trunk 

-corrtab- also allows surpressing correlations based on their p-values when using the option -print(#)-:

corrtab price-trunk, print (.10)

Or, based on the absolute value of the correlation coefficient when specifying -above(#)-:

corrtab price-trunk, above(.4)

Correlations can also be surpressed based on the position in the list of variables. The option -var(#)- only lists the first # variables in the columns, later variables are only included in the rows of the table:

corrtab price-trunk, var(3)

Correlations can be sorted by size for a single variable specified in -vsort()-:

corrtab price-trunk, vsort(price)

In order to use variable labels, some complicated additional commands need to be used; and the default for treating missing values is also pairwise deletion.
In sum, there seems to be little reason to ever use -corrtab-.

 

-makematrix-

-makematrix- by Nick Cox is a very flexible tool for all sort of applications, one of them being correlation tables. He explains it a bit in Cox (2003), but it seems to require too much tweaking to get readymade publication-style correlation tables in a short amount of time. Furthermore, there's a bug in the command: The option -listwise- will yield a pairwise correlation matrix instead of a correlation based on listwise deletion.

In sum, the commands reviewed here all leave a lot to be desired; -mkcorr- is the command that seems to be the most useful one to me.

References

Cox, Nicholas J. 2003. "Speaking Stata. Problems with Tables, Part II." Stata Journal 3(4):420-439.