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
Showing posts with label esttab. Show all posts
Showing posts with label esttab. Show all posts
Mar 24, 2023
Clausen (1998): Applied Correspondence Analysis
Nov 24, 2021
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")
May 7, 2020
Centering predictor variables in OLS regression
clear clear set seed 1 set obs 1000 // Generate variables generate female = (runiform() > .48) generate health = 1 - female * rnormal(0.5, 1) + 0.1 * rnormal(3, 1) // Recode predictor variable center female recode female (0 = -1), gen(e_female) eststo clear eststo: regress health female eststo: regress health c_female eststo: regress health e_female esttab, cells(b(fmt(2))) rename(c_female female e_female female) /// mtitle("Dummy" "Centering" "Effect") nonumber /// coeflabel(female "Female (ref. male)" _cons "Intercept") /// varwidth(18) collabel("") stats(N, fmt(%9.0gc)) /// title(Comparing different types of centering)
Nov 25, 2019
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)
Oct 18, 2019
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
Oct 11, 2019
Indent value labels of categorical variables for -esttab-
// Open Allbus cumulation
use hs01 iscd975 age sex ///
using "ZA4586_v1-0-0.dta", clear
// Self-rated health
recode hs01 (4 5 = 1 "Poor health") ///
(3 2 1 = 0 "Good health") ///
(-1 -11 -9 = .) ///
, gen(poorhealth)
label var poorhealth "Poor self-rated health"
// Education
recode iscd975 (5 = 5 "Tertiary") ///
(4 = 4 "Post-secondary") ///
(3 = 3 "Upper secondary") ///
(2 = 2 "Lower secondary") ///
(1 = 1 "Basic") (-32 = .) ///
, gen(isced)
label var isced "Education"
// Age
recode age (-32 = .)
generate age10 = age/10
label var age10 "Age / 10"
// Sex
recode sex (1 = 0 "Male") ///
(2 = 1 "Female") ///
, gen(female)
label var female "Female sex"
// Indent value labels
labvalch3 female isced, prefix("\dumm ")
// Add this to Latex code
// % Indent for tables
// \newcommand*{\dumm}{\hspace*{0.5cm}}%
eststo clear
eststo: regress poorhealth ib3.isced age10 i.female, robust
esttab using table.tex, b(2) se(2) booktabs replace ///
label nobaselevels nonumber ///
modelwidth(25) varwidth(33) ///
stats(N, fmt(%8.0gc)) ///
refcat(1.isced "Education (ref. Upper secondary)" ///
1.female "Sex (ref. Male)", nol) ///
title("Poor self-rated health regressed on education, age, and sex, linear probability model")
// Alternative approach (that also works for .rtf tables)
labvalch3 female isced, subst("\dumm " "")
*ssc install elabel
elabel define female isced (= #) (= " " + @) , modify
// Note that the leading sapce don't seem to show up in a frequency table with -fre-
esttab using table.rtf, b(2) se(2) replace ///
label nobaselevels nonumber ///
modelwidth(25) varwidth(33) ///
stats(N, fmt(%8.0gc)) ///
refcat(1.isced "Education (ref. Upper secondary)" ///
1.female "Sex (ref. Male)", nol) ///
title("Poor self-rated health regressed on education, age, and sex, linear probability model")
Sep 4, 2018
Making tables for logit models using -esttab-
// Open Allbus 2016 use "ZA5250_v2-0-0.dta", clear // Prepare variables recode hs01 (4 5 = 1 "Poor health") (3 2 1 = 0 "Good health") (-9 = .), gen(poorhealth) label var poorhealth "Poor self-rated health" recode sex (2 = 1 " Female") (1 = 0 " Male"), gen(female) label var female "Female sex" recode age (-32 = .) label var age "Age" recode isced97 (1 2 = 0 " Low") (3 4 = 1 " Medium") (5 6 = 2 " High") (-32 = .), gen(education) label var education "Education" // Fit model 1 eststo clear eststo: qui logit poorhealth i.female i.education age estadd expb qui sum poorhealth estadd scalar avg = r(mean) * 100 // Long table
esttab, cells(b(star fmt(2) label("B")) ///
se(par fmt(2) label("(SE B)")) ///
expb(par([ ]) label("[OR]"))) ///
stats(N avg chi2 df_m, fmt(%8.0gc %8.1f %8.1gc 0) ///
label("Observations" ///
"% poor health" ///
"Chi-squared" ///
"df")) ///
label varwidth(30) modelwidth(25) nonumber nomtitle varlabel(_cons "Intercept") ///
eqlabel(" ") ///
nobaselevel ///
refcat(1.female "Sex (ref. male)" 1.education "Education (ref. low)", nol) ///
addnote("* p<0.05, ** p<0.01, *** p<0.001") ///
title("Poor self-rated health regressed on sex, education, and age. Logistic model")
// Fit models 2
eststo clear
eststo: qui logit poorhealth i.education age if female == 0
qui estadd expb
qui sum poorhealth if female == 0
qui estadd scalar avg = r(mean) * 100
eststo: qui logit poorhealth age i.education if female == 1
qui estadd expb
qui sum poorhealth if female == 1
qui estadd scalar avg = r(mean) * 100
// Wide table
esttab, cell("b(fmt(2) label(B)) se(fmt(2) label(SE B)) expb(fmt(2) label(OR) star)") ///
stats(N avg chi2 df_m, fmt(%8.0gc %8.1f %8.1gc 0) ///
label("Observations" ///
"% poor health" ///
"Chi-squared" ///
"df")) ///
label varwidth(30) modelwidth(8) nonumber mtitle("Men" "Women") varlabel(_cons "Intercept") ///
eqlabel(" ") nobaselevels ///
refcat(1.education "Education (ref. low)", nol) ///
addnote("* p<0.05, ** p<0.01, *** p<0.001") ///
title("Poor self-rated health regressed on education and age, stratified by sex. Logistic model")
Aug 2, 2018
Meta-regression
// Open ISSP 2005 use COUNTRY WRKHRS V51 using "ZA4350_v2-0-0.dta", clear // Prepare variables rename COUNTRY country // Job satisfaction generate jobsatis = 7 - V51 rename WRKHRS workhrs // Average working hours per country bys country: egen avghrs = mean(workhrs) // Center working hours center workhrs, inplace // Fit models eststo clear eststo: regress jobsatis workhrs avghrs // OLS regression eststo: mixed jobsatis workhrs avghrs || country: // RE regression qui matrix foo = e(N_g) // Number of individuals qui estadd scalar nc = foo[1,1] // Number of individuals // Prepare data for meta-regression preserve // Generate country-level data set of average working hours egen pickone = tag(country) keep if pickone keep country avghrs tempfile temp save `temp', replace restore // Average job satisfaction per country, controlling for individual working hours statsby _b[_cons] _se[_cons] e(N), clear by(country): regress jobsatis workhrs // Merge with average working hours merge 1:1 country using `temp' eststo: regress _stat_1 avghrs // OLS regression two-step eststo: vwls _stat_1 avghrs, sd(_stat_2) // FE meta-regression eststo: metareg _stat_1 avghrs, wsse(_stat_2) // RE meta-regressionesttab, keep(main:) mtitles("OLS" "Multilevel" "OLS two-step" "FE meta-regression" "RE meta-regression") b(3) se(3) /// coeflabel(workhrs "Working hours" /// avghrs "Average working hours" /// _cons "Intercept") /// stats(nc N, labels("No. countries") fmt(%12.0gc %12.0gc)) varwidth(25) modelwidth(20) eqlabels("", none)
Jun 10, 2018
Direct and indirect effects in OLS regression
sysuse auto, clear eststo clear // Bivariate relationship: eststo: regress price mpg // Accounting for confounding variable: eststo: regress price mpg displacement // Accounting for confounding variable in three steps: // 1 quietly regress mpg displacement predict u, residuals // 2 quietly regress price displacement predict v, residuals // 3 eststo: regress v u // Table: esttab, rename(u mpg) se nomtitle label varwidth(30)
Mar 11, 2018
Allison's (2014) book on event history and survival analysis
This replicates the analyses in Allison (2014). There is another approach here, but I like mine better.
// Table 2.1
use https://statisticalhorizons.com/wp-content/uploads/rank.dta, clear
ltable dur promo, failure hazard noadjust
// Table 2.2
use https://statisticalhorizons.com/wp-content/uploads/rank.dta, clear
generate id = _n
reshape long art cit, i(id) j(year)
drop if year > dur // Remove empty observation created during reshape
replace promo = 0 if year < dur // Create time-varying failure variable
generate jobpres = prest1 // Create-time varying prestige varianle
replace jobpres = prest2 if year >= jobtime
eststo clear
eststo: logit promo undgrad phdmed phdprest jobpres art cit
estadd expb
eststo: logit promo undgrad phdmed phdprest jobpres art cit year c.year#c.year
estadd expb
esttab, cell("b(fmt(3)) z(fmt(2) star) expb(fmt(2))") varwidth(17) label nonumber ///
mtitle("Model 1" "Model 2") interaction( X ) ///
stats(ll N, label("Log-likelihood") fmt(2 %9.0gc)) ///
title("Logistic Models Predicting the Probability of Promotion") ///
legend
// Likelihood ratio test on p. 13
lrtest est1 est2
// Additional model on p. 13
logit promo undgrad phdmed phdprest jobpres art cit year c.year#c.year c.phdprest#c.year
// Table 2.3
// Scenario A: All censored cases are in fact promoted
use https://statisticalhorizons.com/wp-content/uploads/rank.dta, clear
generate id = _n
reshape long art cit, i(id) j(year)
drop if year > dur // Remove empty observation created during reshape
replace promo = 0 if year < dur // Create time-varying failure variable
replace promo = 1 if year == dur & dur < 10 // All censored cases are promoted
generate jobpres = prest1 // Create-time varying prestige varianle
replace jobpres = prest2 if year >= jobtime
eststo: logistic promo undgrad phdmed phdprest jobpres art cit year c.year#c.year
// Scenario B: All censored cases did not experience event until the end of the observation period
use https://statisticalhorizons.com/wp-content/uploads/rank.dta, clear
generate id = _n
reshape long art cit, i(id) j(year)
drop if year > dur & promo == 1 // Remove empty observation created during reshape
// if they are promoted
replace promo = 0 if year < dur // Create time-varying failure variable
replace art = art[_n-1] if art == . // Carry observations forward
replace cit = cit[_n-1] if cit == . // Carry observations forward
generate jobpres = prest1 // Create-time varying prestige varianle
replace jobpres = prest2 if year >= jobtime
eststo: logistic promo undgrad phdmed phdprest jobpres art cit year c.year#c.year
esttab est2 est3 est4, eform b(2) z wide staraux varwidth(17) label nonumber ///
mtitle("Standard analysis" "Scenario A" "Scenario B") interaction( X ) ///
stats(ll N, label("Log-likelihood") fmt(2 %9.0gc)) ///
title("Extreme Case Scenarios for Informative Censoring") ///
legend modelwidth(18)
// Table 3.1
use https://statisticalhorizons.com/wp-content/uploads/recid.dta, clear
stset week, failure(arrest==1)
eststo clear
eststo: streg fin age race wexp mar paro prio, dist(exponential)
estadd expb
eststo: streg fin age race wexp mar paro prio, dist(weibull)
estadd expb
eststo: streg fin age race wexp mar paro prio, dist(ggamma)
estadd expb
esttab, cell("b(fmt(3)) z(fmt(2) star) expb(fmt(3))") varwidth(17) label nonumber ///
mtitle("Exponential" "Weibull" "Gamma") ///
stats(ll N, label("Log-likelihood") fmt(2 %9.0gc)) ///
title("Estimates for Three Models of Recidivism") ///
legend keep(_t:*) eqlabels("", none) /// // Removes equation label
coeflabel(fin "Financial aid" age "Age at release" race "Black" ///
wexp "Work experience" mar "Married" paro "Paroled" ///
prio "Prior convictions")
// Table 3.2
use https://statisticalhorizons.com/wp-content/uploads/recid.dta, clear
stset week, failure(arrest==1)
eststo clear
qui eststo: streg fin age race wexp mar paro prio, dist(ggamma)
qui estadd scalar dev = -2*e(ll)
qui eststo: streg fin age race wexp mar paro prio, dist(lognormal)
qui estadd scalar dev = -2*e(ll)
qui eststo: streg fin age race wexp mar paro prio, dist(llogistic)
qui estadd scalar dev = -2*e(ll)
qui eststo: streg fin age race wexp mar paro prio, dist(weibull)
qui estadd scalar dev = -2*e(ll)
qui eststo: streg fin age race wexp mar paro prio, dist(gompertz)
qui estadd scalar dev = -2*e(ll)
qui eststo: streg fin age race wexp mar paro prio, dist(exponential)
qui estadd scalar dev = -2*e(ll)
qui esttab, cells(none) scalars("dev Deviance" "aic AIC" "bic BIC") sfmt(3) nomtitles noobs
esttab r(stats, transpose), coeflabel(est1 "Gamma" est2 "Log-normal" ///
est3 "Log-logistic" est4 "Weibull" ///
est5 "Gompertz" est6 "Exponential") ///
title("Goodness of Fit for Recidivism Models") ///
nomtitle collabels("Deviance" "AIC" "BIC")
// Figure 3.1: Hazard Function for Weibull Regression Model
qui streg fin age race wexp mar paro prio, dist(weibull)
stcurve, hazard xtitle("Weeks since release")
// Figure 3.2: Hazard Function for Log-Logistic Regression Model
qui streg fin age race wexp mar paro prio, dist(llogistic)
stcurve, hazard xtitle("Weeks since release")
// Table 4.1
use https://statisticalhorizons.com/wp-content/uploads/recid.dta, clear
stset week, failure(arrest==1)
eststo clear
eststo: stcox fin age race wexp mar paro prio
estadd expb
generate id = _n
reshape long work, i(id) j(stop)
generate start = stop - 1
drop if stop > week // Drop empty observations
replace arrest = 0 if week != stop // Create time-varying failure variable
stset stop, failure(arrest == 1) id(id) origin(start)
eststo: stcox fin age race wexp mar paro prio work
estadd expb
generate worklag = work[_n-1] if start > 0
eststo: stcox fin age race wexp mar paro prio worklag
estadd expb
esttab, cell("b(fmt(3)) z(fmt(2) star) expb(fmt(3))") varwidth(18) label nonumber ///
mtitle("Basic" "Time-varying X" "Lagged X") modelwidth(15) ///
stats(ll N, label("Log-likelihood") fmt(2 %9.0gc)) ///
title("Cox Regression Estimates for Recidivism Data") legend ///
coeflabel(fin "Financial aid" age "Age at release" race "Black" ///
wexp "Work experience" mar "Married" paro "Paroled" ///
prio "Prior convictions" work "Employment") rename(worklag work)
// Table 4.2
list id stop start arrest work fin age if inlist(id, 339, 417), noobs sepby(id)
// Table 4.3
estimates restore est2
estat phtest, detail
// Table 4.4
eststo clear
eststo: stcox fin age race wexp mar paro prio work, tvc(age wexp)
estadd expb
esttab, cell("b(fmt(3)) z(fmt(2) star) expb(fmt(3))") varwidth(22) label nonumber ///
mtitle("Interactions with time") modelwidth(22) ///
stats(ll N, label("Log-likelihood") fmt(2 %9.0gc)) ///
title("Cox Regression Estimates with Time Interactions") legend ///
eqlabels("Main effects" "Interactions with time") ///
coeflabel(fin "Financial aid" age "Age at release" race "Black" ///
wexp "Work experience" mar "Married" paro "Paroled" ///
prio "Prior convictions" work "Employment")
eststo clear
eststo: stcox fin age race wexp mar paro prio work, strata(wexp)
estadd expb
esttab, cell("b(fmt(3)) z(fmt(2) star) expb(fmt(3))") varwidth(22) label nonumber ///
mtitle("Stratification") modelwidth(22) ///
stats(ll N, label("Log-likelihood") fmt(2 %9.0gc)) ///
title("Cox Regression Estimates with Stratification") legend ///
coeflabel(fin "Financial aid" age "Age at release" race "Black" ///
wexp "Work experience" mar "Married" paro "Paroled" ///
prio "Prior convictions" work "Employment") dropped(--)
// Table 4.5
qui stcox fin age race wexp mar paro prio work, tvc(age wexp)
foreach x of numlist 0 10 20 30 40 50 {
local b1 = _b[age] + `x' * _b[tvc:age]
local expb1 = exp(`b1')
local b2 = _b[wexp] + `x' * _b[tvc:wexp]
local expb2 = exp(`b2')
matrix results = (`x' , `b1' , `expb1' , `b2', `expb2')
if `x' == 0 matrix table45 = results
else matrix table45 = (table45\results)
}
esttab matrix(table45, fmt(0 3 3 3 3)), nomtitle ///
collabel("Weeks" "b" "Exp(b)" "b" "Exp(b)") varwidth(0) ///
title("Effects of Age and Work Experience at Different Times")
// Table 4.6
use https://statisticalhorizons.com/wp-content/uploads/recid.dta, clear
stset week, failure(arrest==1)
stcox fin age race wexp mar paro prio
stcurve, survival at(fin = 1 age = 21 race = 1 wexp = 1 mar = 0 paro = 1 prio = 4) ///
outfile(surv, replace)
use surv, clear
egen pickone = tag(_t)
list _t surv1 if inlist(_t, 0, 5, 10, 15, 20, 25, 30, 35, 40, 45, 50, 52) & pickone, noobs sep(0)
// Table 5.1
use https://statisticalhorizons.com/wp-content/uploads/tarp.dta, clear
eststo clear
stset arrstday, failure(type == 1 2)
eststo: stcox fin age white male married paro numprop crimprop numarst edcomb
estadd expb
stset arrstday, failure(type == 1)
eststo: stcox fin age white male married paro numprop crimprop numarst edcomb
estadd expb
stset arrstday, failure(type == 2)
eststo: stcox fin age white male married paro numprop crimprop numarst edcomb
estadd expb
esttab, cell("b(fmt(3)) z(fmt(2) star) expb(fmt(3))") varwidth(22) label nonumber ///
mtitle("All arrests" "Property arrests" "Non-property arrests") modelwidth(18) ///
stats(ll N, label("Log-likelihood") fmt(2 %9.0gc)) ///
title("Estimates of Proportional Hazards Models for Different Arrest Types") legend ///
coeflabel(fin "Financial aid" age "Age at release" white "White" ///
male "Male" married "Married" paro "Paroled" ///
numprop "No. of property convictions" ///
crimprop "Imprisoned for property crime" ///
numarst "No. of arrests" ///
edcomb "Education")
// Figure 5.1
stset arrstday, failure(type == 1)
stcompet cumin = ci, compet1(2)
sort _t
twoway (line cumin _t if type == 1) ///
(line cumin _t if type == 2) ///
(scatteri .2 375 "Property", msymbol(none) mlabpos(0)) ///
(scatteri .15 370 "Non-property", msymbol(none) mlabpos(0)) ///
, legend(off) ///
xtitle("Days since release") ytitle("Probability of arrest") ///
title("Cumulative incidence")
// Table 5.2
eststo clear
stset arrstday, failure(type == 1 2)
eststo: stcox fin age white male married paro numprop crimprop numarst edcomb
estadd expb
stset arrstday, failure(type == 2)
eststo: stcrreg fin age white male married paro numprop crimprop numarst edcomb, compete(type == 1)
estadd expb
stset arrstday, failure(type==1)
eststo: stcrreg fin age white male married paro numprop crimprop numarst edcomb, compete(type==2)
estadd expb
esttab, cell("b(fmt(3)) z(fmt(2) star) expb(fmt(3))") varwidth(22) label nonumber ///
mtitle("All arrests" "Property arrests" "Non-property arrests") modelwidth(18) ///
stats(ll N, label("Log-likelihood") fmt(2 %9.0gc)) ///
title("Estimates of Subdistribution Hazards Models for Different Arrest Types") legend ///
coeflabel(fin "Financial aid" age "Age at release" white "White" ///
male "Male" married "Married" paro "Paroled" ///
numprop "No. of property convictions" ///
crimprop "Imprisoned for property crime" ///
numarst "No. of arrests" ///
edcomb "Education")
// Figure 5.2
qui stcrreg fin age white male married paro numprop crimprop numarst edcomb, compete(type==2)
stcurve, cif xtitle("Days since release") ///
title("Cumulative incidence of non-property arrests") ///
ylab(, format(%6.2f))
// Table 6.1
use https://statisticalhorizons.com/wp-content/uploads/tarp.dta, clear
eststo clear
estpost tabulate arrstcount
esttab, cell(b) nonumber collabel("Number of persons", lhs("Number of arrests")) nomtitle noobs ///
modelwidth(20) varwidth(20) varlabels(, blist(Total "{hline @width}{break}")) ///
title(Frequency Distribution for Number of Arrests)
// Table 6.2
eststo clear
eststo: nbreg arrstcount fin age white male married paro numprop crimprop numarst edcomb
estadd expb
use https://statisticalhorizons.com/wp-content/uploads/arrests.dta, clear
stset length, failure(arrind == 1)
eststo: stcox fin age white male married paro numprop crimprop numarst edcomb
estadd expb
stcox fin age white male married paro numprop crimprop numarst edcomb, cluster(id)
set matsize 942
eststo: stcox fin age white male married paro numprop crimprop numarst edcomb, shared(id)
estadd expb
esttab, cell("b(fmt(3)) z(fmt(2) star) expb(fmt(3))") varwidth(22) label nonumber ///
mtitle("Negstive binomial count model" "Cox regression, gap time" "Cox regression, shared frailty") modelwidth(18) ///
stats(ll N, label("Log-likelihood") fmt(2 %9.0gc)) ///
title("Regression Models for Repeated Arrests") legend ///
coeflabel(fin "Financial aid" age "Age at release" white "White" ///
male "Male" married "Married" paro "Paroled" ///
numprop "No. of property convictions" ///
crimprop "Imprisoned for property crime" ///
numarst "No. of arrests" ///
edcomb "Education") keep(main:*) drop(_cons)
// Table 6.3
use https://statisticalhorizons.com/wp-content/uploads/arrests.dta, clear
stset length, failure(arrind == 1)
eststo clear
eststo: streg spellnum fin age white male married paro numprop crimprop numarst edcomb, cluster(id) dist(weibull)
estadd expb
eststo: streg spellnum fin age white male married paro numprop crimprop numarst edcomb, shared(id) dist(weibull)
estadd expb
stset end, failure(arrind==1) origin(begin)
eststo: stcox fin age white male married paro numprop crimprop numarst edcomb, cluster(id)
estadd expb
esttab, cell("b(fmt(3)) z(fmt(2) star) expb(fmt(3))") varwidth(22) label nonumber ///
mtitle("Weibull, robust z" "Weibull, shared frailty" "Cox regression, origin times") modelwidth(18) ///
stats(ll N, label("Log-likelihood") fmt(2 %9.0gc)) ///
title("Regression Models for Repeated Arrests") legend ///
coeflabel(fin "Financial aid" age "Age at release" white "White" ///
male "Male" married "Married" paro "Paroled" ///
numprop "No. of property convictions" ///
crimprop "Imprisoned for property crime" ///
numarst "No. of arrests" ///
edcomb "Education") keep(main:)
// Analyses on p. 74
stcox fin age white male married paro numprop crimprop numarst edcomb, cluster(id) tvc(numarst) texp(_t/30.4)
Reference
Allison, Paul D. 2014. Event History and Survival Analysis, 2nd ed. Sage. doi: 10.4135/9781452270029
Feb 15, 2018
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
Feb 14, 2018
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
Feb 13, 2018
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
Feb 12, 2018
Interactions in multiple imputation
This replicates the analyses for Table 6.1 for Allison (2002).
// Table 6.1
// Method 1
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
eststo clear
eststo: mi estimate, post: regress gradrat lenroll i.private##c.csat stufac rmbrd
// Method 2
use "https://statisticalhorizons.com/wp-content/uploads/college.dta", clear
mi set mlong
mi register imputed gradrat csat lenroll stufac rmbrd act
mi impute mvn gradrat csat lenroll stufac rmbrd act, ///
add(5) burnin(500) burnbetween(200) emlog emoutput ///
by(private)
eststo: mi estimate, post: regress gradrat lenroll i.private##c.csat stufac rmbrd
// Method 3
use "https://statisticalhorizons.com/wp-content/uploads/college.dta", clear
generate privateXcsat = private * csat
mi set mlong
mi register imputed gradrat csat private lenroll stufac rmbrd act privateXcsat
mi impute mvn gradrat csat private lenroll stufac rmbrd act privateXcsat, ///
add(5) burnin(500) burnbetween(200) emlog emoutput
eststo: mi estimate, post: regress gradrat lenroll i.private csat privateXcsat stufac rmbrd
esttab, not p wide nostar noobs varlabel(_cons "Intercept") ///
order(_cons csat lenroll stufac 1.private rmbrd) ///
rename(privateXcsat 1.private#c.csat) varwidth(25) nobaselevels ///
title(Regression with interaction terms--three methods) ///
mtitle("Method 1" "Method 2" "Method 3") nonumbers
Reference
Allison, Paul D. 2002. Missing Data. Sage. doi: 10.4135/9781412985079
Feb 9, 2018
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
Jan 11, 2018
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", cleareststo 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
Jan 10, 2018
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.1eststo 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
Jan 9, 2018
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
Aug 22, 2017
Multiple imputation in Stata
clear
// Generate data
set obs 10000
generate e = rnormal()
drawnorm x1 x2, ///
corr(1, .4, 1) cstorage(lower)
generate y = 1 + 1 * x1 + 1 * x2 + 25 * e
drop e
// Model with data-generating process
eststo clear
eststo: regress y x1 x2
// Drop cases conditional on other covariate
replace x1 = . if rnormal() > .1 & x2 > -1
// Model with missing data
eststo: regress y x1 x2
// Impute data
mi set mlong // Declare data format
mi register imputed x1 // Declare variable to be imputed
mi impute chained (regress) x1 = y x2, add(10) // Declare imputation model
// Model with imputed data
eststo: mi estimate, post: regress y x1 x2 // post option is important for -esttab-
// Table
esttab, r2 se mtitles("DGP" "Missing data" "Imputed")
Jul 20, 2017
IV regression using the -sem- command
// Read in data from Angrist and Krueger (1991) // https://economics.mit.edu/faculty/angrist/data1/data/angkru1991 infile lwklywge educ yob qob pob using asciiqob.txt, clear // Generate dummy variables as SEM command does not take factor variables qui tabulate qob, gen(qobx) qui tabulate yob, gen(yobx) qui tabulate pob, gen(pobx) drop qobx1 yobx1 pobx1 // Get rid of reference category eststo clear // Model 2 of Table 4.1.1 of Mostly Harmless Econometrics eststo: regress lwklywge educ yobx* pobx*, robust // Model 6 of Table 4.1.1 of Mostly Harmless Econometrics eststo: ivregress 2sls lwklywge yobx* pobx* (educ = qobx*), robust // Model 6 of Table 4.1.1 using the sem command eststo: sem (lwklywge <- yobx* pobx* educ) (educ <- pobx* qobx*), cov(e.lwklywge*e.educ)esttab, b(3) se(3) nostar drop(_cons educ:) /// indicate("9 year-of-birth dummies = yobx*" /// "50 state-of-birth dummies = pobx*") /// title("OLS and 2SLS estimates of the economic returns to schooling") /// coeflabel(educ "Years of education") /// mtitles("OLS" "-ivregress-" "-sem-") nonumbers varwidth(30) /// eqlabels("", none) // Removes equation label
Subscribe to:
Posts (Atom)