Comparing margins from OLS and IV models

webuse hsng2, clear

sum hsngval // Get avg. and SD

// OLS
regress rent pcturban hsngval faminc i.region 
margins, at(hsngval = (32714 48484 64254)) saving(ols, replace)

// IV
ivregress 2sls rent pcturban (hsngval = faminc i.region)
margins, at(hsngval = (32714 48484 64254)) saving(iv, replace)

// Plot
combomarginsplot ols iv, labels("OLS" "IV") ///
                         title("{bf:A}", span bexpand justification(left)) ///
                         xtitle(Median housing value) ///
                         ytitle(Median gross rent) ///
                         xlabel(48484 "Avg." 64254 "+1 SD" 32714 "-1 SD") ///
                         recastci(rarea) ciopt(color(%30) lcolor(white)) ///
                         legend(pos(11) ring(0)) ///
                         scheme(lean1) name(ols_iv_predictions, replace)

Mummolo & Peterson (2018): Fixed-effects regression

set scheme lean1

webuse nlswork, clear
xtset idcode year

// Model of interest:
xtreg ln_w ttl_exp age c.age#c.age tenure c.tenure#c.tenure not_smsa south i.year, fe
local b : display %4.2f _b[ttl_exp]

// Isolate relevant variation in the treatment:
reghdfe ttl_exp, absorb(idcode year) residuals(relevantvariation)

// Identify a plausible counterfactual shift in X given the data:
center ttl_exp

sum relevantvariation
local sd1 : display %4.2f r(sd)

sum ttl_exp
local sd2 : display %4.2f r(sd)

twoway (kdensity c_ttl_exp) ///
       (kdensity relevantvariation), ///
       ytitle("Density") xtitle(Treatment (centered around 0)) ///
       title("{bf:A} Treatment distribution", span bexpand justification(left)) ///
       legend(order(1 "Original distribution, {it:SD} = `sd2'" ///
                    2 "Distribution after FE, {it:SD} = `sd1'" ) ///
                    pos(2) ring(0)) ///
       name(figure1A, replace)

      // Calculate range within individuals
bysort idcode: egen max = max(ttl_exp)
bysort idcode: egen min = min(ttl_exp)
generate range = max - min 
drop max min

qui sum range, detail
local mean : display %4.2f r(mean)
local p95 :  display %4.2f r(p95)

twoway (histogram range, freq) ///
       (scatteri 2000 `mean' (3) "Mean = `mean'", msymbol(i)) ///
       (scatteri 2000 `p95' (3) "95th percentile = `p95'" , msymbol(i)), ///
       title("{bf:B} Within-individual range of treatment", span bexpand justification(left)) ///
       ylabel(, format(%7.0gc)) ytitle(Frequency) ///
       xline(`mean') xline(`p95') ///
       xtitle(Within-individual range) legend(off) ///
       name(figure1B, replace)

graph combine figure1A figure1B, col(2)

// Calculate % of cases where treatment does not vary
egen pickone = tag(idcode)

fre range if pickone
  // In 12% of cases the treatment does not vary

// Multiply the estimated coefficient of interest by the revised standard deviation

display "Estimated coefficient of interest = `b'"
display "Within-standard deviation = `sd1'"
display "Product = " _skip(1) `sd1' * `b'

Reference Mummolo, Jonathan, and Erik Peterson. 2018. "Improving the Interpretation of Fixed Effects Regression Results." Political Science Research and Methods 6(4):829-835. doi: 10.1017/psrm.2017.44

Clausen (1998): Applied Correspondence Analysis

clear all

// Table 1.1
  // Create Table -- the -camat- command allows working with tabular data right 
  //                 away
input leisure class weight
 1 1 301
 1 2 497
 1 3 208
 1 4  50
 1 5 254
 1 6 187
 2 1 261
 2 2 550
 2 3 250
 2 4  27
 2 5 339
 2 6 157
 3 1 361 
 3 2 534
 3 3 204
 3 4  59
 3 5 324
 3 6 216
 4 1 463 
 4 2 766
 4 3 334
 4 4  72
 4 5 350
 4 6 601
 5 1  89
 5 2 350
 5 3 195
 5 4  12
 5 5 143
 5 6 167
 6 1  23 
 6 2 182
 6 3 124
 6 4  10
 6 5  60
 6 6 110
 7 1 117
 7 2 298
 7 3 145
 7 4  11
 7 5 184
 7 6  56
 8 1 104
 8 2 379
 8 3 219
 8 4  21
 8 5 152
 8 6 213
 9 1 130
 9 2 352
 9 3 153
 9 4  17
 9 5 272
 9 6 264
10 1 168
10 2 370
10 3 187
10 4  51
10 5 162
10 6 424
end

label define leisure 1 "Sport events" ///
                     2 "Cinema" ///
                     3 "Dance/disco" ///
                     4 "Cafe/restaurant" ///
                     5 "Theater" ///
                     6 "Classical concert" ///
					 7 "Pop concert" ///
					 8 "Art exhibition" ///
					 9 "Library" ///
					10 "Church service"
					 
label define class 1 "Manual" ///
                   2 "Low nonmanual" ///
                   3 "High nonmanual" ///
                   4 "Farmer" ///
                   5 "Student" ///
                   6 "Retired" 

label values leisure leisure
label values class class

label var leisure "Leisure activities"
label var class "Occupational class"

expand weight

estpost tabulate leisure class
esttab, unstack nonumber varwidth(20) compress nonote noobs

// Figure 1.1
ca leisure class, plot

// Make a nicer version of Figure 1.1 
capture frame drop biplotc
frame create biplotc
frame biplotc: matrix dim = e(TC)
frame biplotc: svmat2 dim, rname(varname) name(col)
frame biplotc: generate variable = "Social class"

capture frame drop biplotr
frame create biplotr
frame biplotr: matrix dim = e(TR)
frame biplotr: svmat2 dim, rname(varname) name(col)
frame biplotr: generate variable = "Leisure activities"

frame biplotc: save deleteme, replace
frame biplotr: append using deleteme
frame biplotc: erase deleteme.dta
frame drop biplotc
capture frame drop biplot
frame rename biplotr biplot

frame biplot: drop dim3-dim5
frame biplot: replace varname = subinstr(varname, "_", " ", .)
frame biplot: twoway (scatter dim2 dim1 if variable == "Leisure activities", mlabel(varname)) ///
                     (scatter dim2 dim1 if variable == "Social class", mlabel(varname)) ///
                     , legend(order(1 "Leisure activities" 2 "Social classes")) ///
                       xtitle(Dimension 1: Young versus old) ///
                       ytitle(Dimension 2: Art versus light entertainment) ///
                       name(figure11, replace) xscale(range(-1 1)) yscale(range(-1 1))

// Table 2.1
clear all

input region crime weight
1 1  395
1 2 2456
1 3 1758
2 1  147
2 2  153
2 3  916
3 1  694
3 2  327
3 3 1347
end 

label define region 1 "Oslo area" 2 "Mid-Norway" 3 "North-Norway"
label define crime 1 "Burglary" 2 "Fraud" 3 "Vandalism"

label val region region 
label val crime crime

label var region "Region"
label var region "Type of crime"
expand weight

// Table 2.2
table region crime, statistic(proportion, across(crime)) nformat(%6.3f)
table region crime, statistic(proportion, across(region)) nformat(%6.3f)

// Figure 2.2
ca region crime, plot

// Make a nicer version of Figure 2.2 
capture frame drop biplotc
frame create biplotc
frame biplotc: matrix dim = e(TC)
frame biplotc: svmat2 dim, rname(varname) name(col)
frame biplotc: generate variable = "Region"

capture frame drop biplotr
frame create biplotr
frame biplotr: matrix dim = e(TR)
frame biplotr: svmat2 dim, rname(varname) name(col)
frame biplotr: generate variable = "Type of crime"

frame biplotc: save deleteme, replace
frame biplotr: append using deleteme
frame biplotc: erase deleteme.dta
frame drop biplotc
capture frame drop biplot
frame rename biplotr biplot

frame biplot: list
frame biplot: replace varname = subinstr(varname, "_", " ", .)
frame biplot: twoway (scatter dim2 dim1 if variable == "Type of crime", mlabel(varname)) ///
                     (scatter dim2 dim1 if variable == "Region", mlabel(varname)) ///
                     , legend(order(1 "Type of crime" 2 "Region")) ///
                       xtitle(Dimension 1) ytitle(Dimension 2) ///
                       name(figure22, replace) xscale(range(-1 1)) yscale(range(-1 1))

// Table 2.4
// Eigenvalues
di "           " _skip(5) "Dim. 1    "       _skip(5) "Dim. 2   "_skip(5)  "Sum"
di "Eigenvalues" _skip(5) (e(Sv)[1,1])^2 _skip(5) (e(Sv)[1,2])^2 _skip(5) e(inertia)

// Figure 2.4
matrix dense = (749, 66 \ 235, 135 \ 283, 185)
matrix colnames dense = "Dense" "Sparse"
matrix rownames dense = oslo midnorway northnorway
matrix norway = (4558, 5129, 10842)
matrix rowname norway = Norway
ca region crime, plot colsupp(dense) rowsupp(norway)


// Make a nicer version of Figure 2.4
capture frame drop biplotc
frame create biplotc
frame biplotc: matrix dim = e(TC)
frame biplotc: svmat2 dim, rname(varname) name(col)
frame biplotc: generate variable = "Region"

capture frame drop biplotr
frame create biplotr
frame biplotr: matrix dim = e(TR)
frame biplotr: svmat2 dim, rname(varname) name(col)
frame biplotr: generate variable = "Type of crime"

capture frame drop biplotsuppc
frame create biplotsuppc
frame biplotsuppc: matrix dim = e(TC_supp)
frame biplotsuppc: svmat2 dim, rname(varname) name(col)
frame biplotsuppc: generate variable = "Population density"
frame biplotsuppc: list

capture frame drop biplotsuppr
frame create biplotsuppr
frame biplotsuppr: matrix dim = e(TR_supp)
frame biplotsuppr: svmat2 dim, rname(varname) name(col)
frame biplotsuppr: generate variable = "Nat'l average"
frame biplotsuppr: list

frame biplotc: save deleteme, replace
frame biplotr: append using deleteme
frame biplotc: erase deleteme.dta
frame biplotsuppc: save deleteme, replace
frame biplotr: append using deleteme
frame biplotsuppc: erase deleteme.dta
frame biplotsuppr: save deleteme, replace
frame biplotr: append using deleteme
frame biplotsuppr: erase deleteme.dta

frame drop biplotc
capture frame drop biplot
frame rename biplotr biplot

frame biplot: list
frame biplot: replace varname = subinstr(varname, "_", " ", .)
frame biplot: twoway (scatter dim2 dim1 if variable == "Type of crime", mlabel(varname)) ///
                     (scatter dim2 dim1 if variable == "Region", mlabel(varname)) ///
                     (scatter dim2 dim1 if variable == "Population density", mlabel(varname)) ///
                     (scatter dim2 dim1 if variable == "Nat'l average", mlabel(varname)) ///
                     (line dim2 dim1 if variable == "Population density")  ///
                     , legend(order(1 "Type of crime" 2 "Region" 3 "Population density")) ///
                       xtitle(Dimension 1) ytitle(Dimension 2) ///
                       name(figure24, replace) xscale(range(-1 1)) yscale(range(-1 1))

// Table 3.1

clear all
input disease age freq
 1 1  12
 1 2  22
 1 3  35
 1 4  68
 1 5 102
 1 6 147
 2 1   7
 2 2  11
 2 3  35
 2 4  45
 2 5  49
 2 6  33
 3 1  44
 3 2  47
 3 3  45
 3 4  42
 3 5  68
 3 6 155
 4 1  12
 4 2   6
 4 3   5
 4 4  38
 4 5 222
 4 6 469
 5 1  63
 5 2  70
 5 3  69
 5 4  74
 5 5  80
 5 6  84
 6 1   0
 6 2   1
 6 3   8
 6 4  14
 6 5  36
 6 6  37
 7 1   9
 7 2   5
 7 3   5
 7 4  15
 7 5  32
 7 6  64
 8 1   8
 8 2   9
 8 3  30
 8 4  28
 8 5  39
 8 6  56
 9 1 103
 9 2 110
 9 3 138
 9 4 124
 9 5  88
 9 6  54
10 1  22
10 2  43
10 3 105
10 4 165
10 5 314
10 6 334
11 1  30
11 2  42
11 3  42
11 4  72
11 5 126
11 6 235
12 1   7
12 2  13
12 3  38
12 4  32
12 5  48
12 6  76
end

label define age 1 "0-6" 2 "7-15" 3 "16-24" 4 "25-44" 5 "45-66" 6 "67+"
label val age age

label define disease  1 "Nervous disorders" ///
                      2 "Nervous system" ///
                      3 "Eye and ear" ///
                      4 "Cardiovascular" /// 
                      5 "Respiratory organ" ///
                      6 "Stomach ulcer" ///
                      7 "Other digestive disease"  ///
                      8 "Urinary/genital system" ///
                      9 "Skin and subcutis" ///
                     10 "Muscuskeletal dis." ///
                     11 "Other diseases" ///
                     12 "Injuries"
label val disease disease

expand freq

// Table 3.2/Figure 3.1
ca disease age, plot

// Make a nicer version of Figure 3.1 
capture frame drop biplotc
frame create biplotc
frame biplotc: matrix dim = e(TC)
frame biplotc: svmat2 dim, rname(varname) name(col)
frame biplotc: generate variable = "Age group"

capture frame drop biplotr
frame create biplotr
frame biplotr: matrix dim = e(TR)
frame biplotr: svmat2 dim, rname(varname) name(col)
frame biplotr: generate variable = "Disease"

frame biplotc: save deleteme, replace
frame biplotr: append using deleteme
frame biplotc: erase deleteme.dta
frame drop biplotc
capture frame drop biplot
frame rename biplotr biplot

frame biplot: list
frame biplot: replace varname = subinstr(varname, "_", " ", .)
frame biplot: twoway (scatter dim2 dim1 if variable == "Disease", mlabel(varname)) ///
                     (scatter dim2 dim1 if variable == "Age group", mlabel(varname)) ///
                     (line dim2 dim1 if variable == "Age group") ///
                     , legend(order(1 "Disease" 2 "Age group")) ///
                       xtitle(Dimension 1) ytitle(Dimension 2) ///
                       name(figure31, replace) xscale(range(-1 1)) yscale(range(-1 1))

					   
// Table 3.4

clear all
input disease age female freq
 1 1 0   8
 1 2 0  25
 1 3 0  23
 1 4 0  46
 1 5 0  66
 1 6 0 106
 2 1 0   4
 2 2 0  13
 2 3 0  30
 2 4 0  29
 2 5 0  29
 2 6 0  31
 3 1 0  38
 3 2 0  52
 3 3 0  41
 3 4 0  48
 3 5 0  72
 3 6 0 156
 4 1 0  10
 4 2 0   7
 4 3 0   5
 4 4 0  39
 4 5 0 226
 4 6 0 399
 5 1 0  70
 5 2 0  77
 5 3 0  72
 5 4 0  77
 5 5 0  89
 5 6 0 100
 6 1 0   0
 6 2 0   0
 6 3 0   8
 6 4 0  22
 6 5 0  43
 6 6 0  43
 7 1 0  12
 7 2 0   3
 7 3 0   3
 7 4 0  16
 7 5 0  21
 7 6 0  60
 8 1 0  10
 8 2 0  11
 8 3 0  11
 8 4 0  12
 8 5 0  22
 8 6 0  72
 9 1 0 105
 9 2 0 106
 9 3 0 119
 9 4 0 106
 9 5 0  72
 9 6 0  55
10 1 0  12
10 2 0  36
10 3 0  95
10 4 0 153
10 5 0 267
10 6 0 266
11 1 0  38
11 2 0  47
11 3 0  41
11 4 0  52
11 5 0 101
11 6 0 192
12 1 0   4
12 2 0  13
12 3 0  55
12 4 0  46
12 5 0  63
12 6 0  77
 1 1 0  17
 1 2 1  19
 1 3 1  45
 1 4 1  89
 1 5 1 135
 1 6 1 179
 2 1 1  10
 2 2 1   9
 2 3 1  41
 2 4 1  60
 2 5 1  68
 2 6 1  34
 3 1 1  50
 3 2 1  44
 3 3 1  48
 3 4 1  36
 3 5 1  64
 3 6 1 155
 4 1 1  15
 4 2 1   5
 4 3 1   5
 4 4 1  38
 4 5 1 218
 4 6 1 523
 5 1 1  56
 5 2 1  64
 5 3 1  67
 5 4 1  70
 5 5 1  72
 5 6 1  72
 6 1 1   0
 6 2 1   3
 6 3 1   8
 6 4 1   5
 6 5 1  29
 6 6 1  33
 7 1 1   6
 7 2 1   6
 7 3 1   6
 7 4 1  15
 7 5 1  42
 7 6 1  68
 8 1 1   6
 8 2 1   8
 8 3 1  48
 8 4 1  45
 8 5 1  55
 8 6 1  43
 9 1 1 100
 9 2 1 113
 9 3 1 156
 9 4 1 141
 9 5 1 104
 9 6 1  54
10 1 1  33
10 2 1  49
10 3 1 115
10 4 1 177
10 5 1 358
10 6 1 387
11 1 1  21
11 2 1  37
11 3 1  44
11 4 1  92
11 5 1 150
11 6 1 268
12 1 1  10
12 2 1  14
12 3 1  23
12 4 1  19
12 5 1  34
12 6 1  74
end

label define age 1 "0-6" 2 "7-15" 3 "16-24" 4 "25-44" 5 "45-66" 6 "67+"
label val age age

label define disease  1 "Nervous disorders" ///
                      2 "Nervous system" ///
                      3 "Eye and ear" ///
                      4 "Cardiovascular" /// 
                      5 "Respiratory organ" ///
                      6 "Stomach ulcer" ///
                      7 "Other digestive disease"  ///
                      8 "Urinary/genital system" ///
                      9 "Skin and subcutis" ///
                     10 "Muscuskeletal dis." ///
                     11 "Other diseases" ///
                     12 "Injuries"
label val disease disease

label define female 1 "Female" 0 "Male"
label val female female
expand freq

// Table 3.5, Table 3.6, Figure 3.2
ca (demo: age female) disease, plot dim(3)
// Figure 3.3
cabiplot, dim(3 2)

// Table 4.1
clear all

input classification subgroup freq
 1 1 139
 1 2  40
 1 3  40
 1 4  41
 2 1 132
 2 2  42
 2 3  37
 2 4  53
 3 1 131
 3 2  21
 3 3  16
 3 4  15
 4 1 124
 4 2  51
 4 3  64
 4 4 124
 5 1 101
 5 2  45
 5 3  49
 5 4  62
 6 1  15
 6 2  79
 6 3   5
 6 4   4
 7 1  20
 7 2  98
 7 3  34
 7 4  29
 8 1  24
 8 2  47
 8 3   1
 8 4   2
 9 1   5
 9 2  42
 9 3  10
 9 4   1
10 1   7
10 2  65
10 3  12
10 4   6
11 1 137 
11 2 114
11 3 106
11 4 159
12 1  61
12 2  67
12 3 115
12 4  62
13 1  95
13 2  44
13 3  83
13 4  86
14 1 143
14 2  83
14 3 121
14 4 149
15 1  57
15 2  97
15 3  92
15 4  98
16 1  76
16 2  32
16 3  56
16 4 195
17 1  75
17 2  49
17 3  43
17 4 194
18 1  63
18 2  45
18 3  38
18 4 171
19 1  48
19 2  46
19 3  18
19 4 143
20 1  49
20 2 113
20 3  46
20 4 105
21 1 111
21 2  11
21 3  76
21 4  95
22 1  21
22 2  38
22 3  23
22 4  34
23 1  24
23 2  36
23 3  33
23 4  60
24 1 115
24 2  50
24 3  66
24 4 106
25 1  64
25 2  90
25 3  66
25 4  69
26 1  24
26 2   8
26 3  15
26 4  72
27 1  71
27 2  26
27 3  40
27 4  58
28 1  57
28 2  32
28 3  29
28 4  82
29 1  72
29 2  53
29 3  55
29 4  65
30 1  30
30 2  37
30 3  37
30 4  52
31 1  25
31 2  44
31 3  18
31 4  28
32 1  86
32 2  15
32 3  66
32 4 149
end

expand freq

label define subgroup 1 "Sick" 2 "Deviant" 3 "Dependent" 4 "Indebted"
label val subgroup subgroup

label define classification ///
 1 "Poor mental health" ///
 2 "Poor general health" ///
 3 "Using sedatives" ///
 4 "National insurance" ///
 5 "Low education" ///
 6 "Alcohol consumption" ///
 7 "Convicted" ///
 8 "Alcohol problems" ///
 9 "Ever used narcotics" ///
10 "Debts due to penalty" ///
11 "Daily cigarette smoking" ///
12 "Long-term client" ///
13 "Trouble daily expenses" ///
14 "Trouble NOK2000" ///
15 "Unemployed" ///
16 "House debt" ///
17 "Owns a dwelling" ///
18 "Owns a car" ///
19 "High income" ///
20 "Men" ///
21 "Women" ///
22 "Age 18-24" ///
23 "Age 25-30" ///
24 "Age 30-50" ///
25 "Unmarried" ///
26 "Married" ///
27 "Divorced" ///
28 "Rural" ///
29 "Urban" ///
30 "City" ///
31 "Child not in household" ///
32 "Child in household"
label val classification classification

ca classification subgroup, plot

// Table 5.1

clear all

input female age alcohol freq
0 0 0  22
0 0 1  78
0 0 2 109
0 0 3 108
0 0 4 132
0 0 5  85
0 1 0  19
0 1 1  84
0 1 2 120
0 1 3  91
0 1 4 203
0 1 5  90
0 2 0  83
0 2 1 130
0 2 2 135
0 2 3 108
0 2 4 160
0 2 5  78
0 3 0  69
0 3 1 126
0 3 2  99
0 3 3  50
0 3 4 127
0 3 5 106
0 4 0  32
0 4 1  62
0 4 2  40
0 4 3  41
0 4 4  95
0 4 5 126
0 5 0  22
0 5 1  19
0 5 2  22
0 5 3  29
0 5 4  62
0 5 5 177
1 0 0  54
1 0 1 134
1 0 2 114
1 0 3 101
1 0 4 104
1 0 5  80
1 1 0  65
1 1 1 127
1 1 2 136
1 1 3  87
1 1 4  81
1 1 5  35
1 2 0 105
1 2 1 150
1 2 2 124
1 2 3  67
1 2 4 113
1 2 5  38
1 3 0 139
1 3 1 149
1 3 2 103
1 3 3  71
1 3 4  88
1 3 5  49
1 4 0  82
1 4 1  63
1 4 2  54
1 4 3  44
1 4 4  84
1 4 5  53
1 5 0  39
1 5 1  40
1 5 2  42
1 5 3  36
1 5 4  72 
1 5 5 100
end
 
label define female 0 "Male" 1 "Female"
label define alcohol 5 "Never"           4 "More seldom" 3 "Once a month" ///
                     2 "2-3 times/month" 1 "Once a week" 0 "Many times a week"
label define age 0 "16-25" 1 "26-35" 2 "36-45" 3 "46-55" 4 "56-66" 5 "67-100"
label val female female
label val alcohol alcohol
label val age age


// Loglinear model
eststo clear
eststo Independent: glm freq i.female  i.age               i.alcohol, fam(pois) link(log)                 // S A C

eststo Base:        glm freq i.female##i.age               i.alcohol, fam(pois) link(log)                 // SA C
eststo M1:          glm freq i.female##i.alcohol    i.female##i.age, fam(pois) link(log)                  // SC SA
eststo M2:          glm freq i.female##i.age        i.age##i.alcohol, fam(pois) link(log)                 // SA AC
eststo M3:          glm freq i.female##i.alcohol    i.age##i.alcohol, fam(pois) link(log)                 // SC AC
eststo M4:          glm freq i.female##i.alcohol    i.age##i.alcohol i.female##i.age, fam(pois) link(log) // SC AC SA
eststo Saturated:   glm freq i.female##i.age##i.alcohol, fam(pois) link(log)                              // SAC

esttab, cells(none) scalars(df deviance deviance_p p) noobs nomtitles nonumber
esttab r(stats, transpose fmt(0 1 1 3)), collabels("df" "L squared" "Pearson chi-squared" "P") ///
                            mtitle("") modelwidth(20) ///
                            labcol2("[S][A][C]" "[SA][C]" "[SC][SA]" "[SA][AC]" "[SC][AC]" "[SC][AC][SA]" "[SAC]")

// Table 5.4
ca (demo: age female) alcohol [fw = freq], plot

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")

Random graphs (148): Dynamite plot

sysuse auto, clear

qui regress mpg i.foreign
local df = e(df_r)
local t = round(((_b[1.foreign]/_se[1.foreign])), .01)
local p = round((2 * ttail(`df', abs(`t'))), .001)
margins, at(foreign = (0 1))
marginsplot, recast(bar) ytitle("Mileage") xtitle("Car type") ///
             plotopts(fcolor(gs14) lcolor(black)) ///
             title("{bf:A} Mileage", justification(left) bexpand span) ///
             addplot(scatteri 28.95 0 28.95 1, recast(line) lwidth(medium) lpattern(solid) ///
                     text(29.5 .5 "{it:t}(`df') = `t', {it:p} = `p'") ///
                     xlabel(0 1) ///
                    || scatteri 28.95 0 28.05 0, recast(line) lwidth(medium) lpattern(solid) ///
                    || scatteri 28.95 1 28.05 1, recast(line) lwidth(medium) lpattern(solid)) ///
             legend(off) ///
             name(figureA, replace)
			 
qui regress weight i.foreign
local df = e(df_r)
local t = round(((_b[1.foreign]/_se[1.foreign])), .01)
local double p = round((2 * ttail(`df', abs(`t'))), .0000001)
margins, at(foreign = (0 1))
marginsplot, recast(bar) ytitle("Weight") xtitle("Car type") ///
             plotopts(fcolor(gs14) lcolor(black)) ///
             title("{bf:B} Weight", justification(left) bexpand span) ///
             addplot(scatteri 3700 0 3700 1, recast(line) lwidth(medium) lpattern(solid) ///
                     text(3760 .5 "{it:t}(`df') = `t', {it:p} = `p'") ///
                     xlabel(0 1) ///
                   || scatteri 3700 0 3600 0, recast(line) lwidth(medium) lpattern(solid) ///
                   || scatteri 3700 1 3600 1, recast(line) lwidth(medium) lpattern(solid)) ///
             legend(off) ///
             name(figureB, replace)
			 
graph combine figureA figureB, col(2) name(figure, replace)

Random graphs (147): Labeling things directly, formatting dates on axes

// Turn date from string to Stata format
gen date2 = date(date, "DMY")
format date2 %td

// Create local that contains the most recent date for labeling lines directly
qui sum date2 
local lastdate = r(max)

// Create list of groups
local groups A B C D

// Create value of outcome for each group 
foreach x of local groups {
	qui sum outcome if group == "`x'"	
	local `x'_count = r(max)  + (log(r(max)) * 10) // This is to offset the label a bit
}

// Figure
twoway (connected outcome date2 if groups == "A", msymbol(x)) ///
       (connected outcome date2 if groups == "B", msymbol(x)) ///
       (connected outcome date2 if groups == "C", lcolor(red) mcolor(red) msymbol(x)) ///
       (connected outcome date2 if groups == "D", msymbol(x)) ///
       (scatteri `A_count' `lastdate' (9) "{it:A}" ///
                 `B_count' `lastdate' (9) "{it:B}" ///
                 `C_count' `lastdate' (9) "{it:C}" ///
                 `D_count' `lastdate' (9) "{it:D}" ///
               , msymbol(i)) ///
      , ytitle("Outcome") ylabel(, format(%9.0gc)) ///
		xtitle("") xlabel(,ang(30) format(%tdMonth_dd)) ///  // This is to format the date
		legend(off) name(figure, replace)

Random graphs (146): Bar graph

clear

// Install rcall command
*github install haghish/rcall, stable

// Clear memory of R session
rcall clear

//#install.packages("essurvey");
rcall: library(essurvey)

// Download data
rcall: download_rounds(c(9), ///
                       ess_email = "REGISTERED E-MAIL ADDRESS HERE", ///
                       output_dir = "data", ///
                       format = 'stata')
					   
// Open data set
local vars cntry gndr polintr dweight
use `vars' using "data/ESS9/ESS9e01_2.dta", clear

// Variables
generate interest = (inlist(polintr, 1, 2))
replace  interest =  . if missing(polintr)
generate female = (gndr == 2)

// Calculate gender gap per country
tempname foo
postfile `foo' str50 cntry men women gap gap_se using deleteme.dta, replace
levelsof cntry, local(cntry)
foreach i of local cntry {
 qui regress interest i.female [pw = dweight] if cntry == `"`i'"'
 qui margins i.female, post
 local   men    = _b[0.female] * 100
 local women    = _b[1.female] * 100
 qui regress interest i.female [pw = dweight] if cntry == `"`i'"'
 qui margins r.female, post
 local   gap    = _b[r1vs0.female] * -100
 local   gap_se = _se[r1vs0.female] * -100
 post `foo' ("`i'") (`men') (`women') (`gap') (`gap_se')
}
postclose `foo'

// Open data set with gender gap per country
use deleteme.dta, clear

kountry cntry, from(iso2c)
egen country = rank(gap), unique
labmask country, value(NAMES_STD)

generate ub = gap + 1.96 * gap_se
generate lb = gap - 1.96 * gap_se

generate y1 = -women
generate y2 = 0 
generate y3 = -women + men

generate x1 = country - 0.2
generate x2 = country + 0.2

mylabels 0(20)60, myscale(-@) local(show)

local spaces = 30 * " "

twoway (scatter y1 country, msymbol(p) yaxis(1 2)) ///
       (rbar y1 y2 x1, blcolor(gs10%40) bfcolor(gs10%40) barwidth(0.4)) ///
       (rbar y1 y3 x2, blcolor(gs5%40) bfcolor(gs5%40) barwidth(0.4)) ///
       (rbar gap y2 x2, blcolor(gs3%40) bfcolor(gs3%40) barwidth(0.4)) ///
       (rspike ub lb x2) ///
     , ylabel(`show',  grid ang(h)) /// 
       ylabel(0(10)40, grid ang(h) axis(2)) /// 
       ytitle("% interested in politics") ///
       ytitle("`spaces' Gender gap", axis(2)) ///
       xlabel(1/19, noticks ang(45) valuelabel labsize(small))  ///
       legend(order(2 "Men" 3 "Women" 4 "Gender gap" 5 "Gender gap 95% CI") ///
              col(4) pos(12) ring(0)) name(gendergap, replace)