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 Textbooks. Show all posts
Showing posts with label Textbooks. 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")
Oct 26, 2019
Specification curve analysis
use "ZA4612_v1-0-1.dta", clear
do "ZA4612_patch_v1-0-1.do" // Some patch from data provider
keep v399-v404 v933 v827 v301 v298
// Outcomes
mvdecode v399-v404, mv(9 = .a )
generate stress = 5 - v399
generate depress = 5 - v400
generate calm = v401 - 1
generate energy = v402 - 1
generate pain = 5 - v403
generate lonely = 5 - v404
// Predictors
// Topbot
mvdecode v933 v827, mv(96 = .a \ 99 = .b)
generate topbot = v933
replace topbot = v827 if missing(topbot)
// Age
mvdecode v301, mv(999 = .a)
rename v301 age
// Female
generate female = (v298 == 2)
drop v399-v404 v933 v827 v298 // Drop unnecessary variables
// Post definition
tempname foo
postfile `foo' str50 spec k1 k2 k3 b se using deleteme.dta, replace
// Loop
forvalues k1 = 1/6 {
if `k1' == 1 local y "stress"
if `k1' == 2 local y "depress"
if `k1' == 3 local y "pain"
if `k1' == 4 local y "calm"
if `k1' == 5 local y "energy"
if `k1' == 6 local y "lonely"
forvalues k2 = 1/4 {
if `k2' == 1 local agecontrol "age"
if `k2' == 2 local agecontrol "c.age##c.age"
if `k2' == 3 local agecontrol "c.age##c.age##c.age"
if `k2' == 4 local agecontrol "c.age##c.age##c.age##c.age"
forvalues k3 = 1/3 {
if `k3' == 1 local ifs " "
if `k3' == 2 local ifs "if female == 1"
if `k3' == 3 local ifs "if female == 0"
local spec regress `y' topbot `agecontrol' female `ifs'
qui `spec'
local b = _b[topbot]
local se = _se[topbot]
post `foo' ("`spec'") (`k1') (`k2') (`k3') (`b') (`se')
}
}
}
postclose `foo'
// Plot specification curve
use deleteme, clear
// Generate ranked analytical choice variable
sort b
generate sk = _n
label var sk "Specification (sorted by coefficient size)"
// Calculate CI's
generate ub = b + 1.96 * se
generate lb = b - 1.96 * se
// Remind yourself what the variables mean
label var k1 "Outcome"
label var k2 "Age control"
label var k3 "Subsample"
// Stack indicators
generate k3c = k3
generate k2c = k2 + 1 + 3 // 3 because K3 has 3 categories, 1 for title
generate k1c = k1 + 1 + 3 + 4 + 1 // K2 has 4 categories
// Calculate some things for size of second axis
qui summarize b
global brange = r(max) - r(min)
global bmin = r(min)
global bmax = r(max)
global from_y = $bmin - (4.5 * $brange)
// Plot
twoway (scatter k1c k2c k3c sk, msymbol(o o o) msize(vsmall vsmall vsmall) ///
yscale(range(1 22)) ///
ylabel( 1 "Full sample" 2 "Females only" 3 "Males only" 4 "{bf:Subsample}" ///
5 "Age" 6 "Age squared" 7 "Age cubed" 8 "Age quartic" 9 "{bf:Age control}" ///
10 "Stress" 11 "Depression" 12 "Pain" ///
13 "Calmness" 14 "Energy" 15 "Loneliness" 16 "{bf:Outcome}", tstyle(notick) axis(1))) ///
(rcap b b sk, yaxis(2) yscale(range($from_y $bmax) axis(2)) ///
ylab(, format(%2.1g) axis(2)) ///
ytitle("{bf:Coefficient size}", axis(2) placement(north)) ///
yline(0, axis(2))) ///
(rspike ub lb sk, yaxis(2)), ///
xlabel(none) ///
legend(off) name(curve, replace)
Reference
Simonsohn, Uri, Joseph P. Simmons, and Leif D. Nelson. 2015. Specification Curve. Descriptive and Inferential Statistics on All Reasonable Specifications. University of Pennsylvania. doi: 10.2139/ssrn.2694998
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
Mar 9, 2018
Mediation analysis
This replicates two illustrations from Iacobucci (2008).
// Section 3.3 (pp. 21-23) // Illustration in section 3.4.1 (pp. 25-26) is identical, just the number of // observations should be 100 // Read in data clear ssd init y x m ssd set observations 50 #delimit ; ssd set corr 1.00 \ 0.45 1.00 \ 0.63 0.55 1.00 ; #delimit cr // SEM approach sem (m <- x) (y <- m x) // Calculation by hand di (_b[m:x] * _b[y:m]) / ((_b[m:x] * _b[y:m]) + _b[y:x]) // Calculation via estat, teffects estat teffects, compact matrix b_indirect = r(indirect) matrix b_total = r(total) scalar indirect = el(b_indirect, 1, 3) scalar total = el(b_total, 1, 3) di "Proportion of total effect mediated by M: " indirect/total // Section 4.2 (pp. 35-38) // Read in data clear ssd init y x m q ssd set observations 50 #delimit ; ssd set corr 1.00 \ 0.45 1.00 \ 0.63 0.55 1.00 \ 0.40 0.40 0.40 1.00 ; #delimit cr // Figure 4.3 // Baseline sem (m <- x) (y <- m x) // (a) sem (m <- x) (y <- m x) (x <- q) // (b) sem (m <- x) (y <- m x) (x -> q) // (d) sem (m <- x) (y <- m x) (m -> q) // (e) sem (m <- x) (y <- m x) (y -> q) // Figure 4.4 // (c) r's = .40 sem (m <- x) (y <- m x) (m <- q) // (e) r's = .40 sem (m <- x) (y <- m x) (y <- q) clear ssd init y x m q ssd set observations 50 #delimit ; ssd set corr 1.00 \ 0.45 1.00 \ 0.63 0.55 1.00 \ 0.70 0.70 0.70 1.00 ; #delimit cr // (c') r's = .70 sem (m <- x) (y <- m x) (m <- q) // (e') r's = .70 sem (m <- x) (y <- m x) (y <- q) //
Reference
Iacobucci, Dawn. 2008. Mediation Analysis. Sage. doi: 10.4135/9781412984966
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 2, 2017
Analyzing natural policy experiments
Hu et al. 2017 simulate data of a natural policy experiment and show how to analyze the data with regression adjustment, propensity score matching, difference-in-differences, and fixed effects regression. (The paper also includes IV, regression discontinuity, and interrupted time series, but does not describe the data created for these analyses.)
clear
set seed 2
// Create data set with experimental conditions
input str4 educ str6 sex str9 treatmentstr str4 health1 number
Low Male Exposed Poor 333
Low Male Exposed Good 917
Low Male Unexposed Poor 1000
Low Male Unexposed Good 2750
Low Female Exposed Poor 500
Low Female Exposed Good 3250
Low Female Unexposed Poor 167
Low Female Unexposed Good 1083
High Male Exposed Poor 83
High Male Exposed Good 542
High Male Unexposed Poor 584
High Male Unexposed Good 3791
High Female Exposed Poor 125
High Female Exposed Good 1750
High Female Unexposed Poor 208
High Female Unexposed Good 2917
end
expand number // Create full number of cases
// Transform strings to numerical variables
generate loeduc = (educ == "Low")
generate female = (sex == "Female")
generate treatment = (treatmentstr == "Exposed")
generate good1 = (health1 == "Good")
label define loeduc 0 "High" 1 "Low"
label val loeduc loeduc
label define female 0 "Male" 1 "Female"
label val female female
label define treatment 0 "Unexposed" 1 "Exposed"
label val treatment treatment
label define health 0 "Poor" 1 "Good"
label val good1 health
// Simulate outcome variable
generate good2 = good1
replace good2 = 1 if good1 == 0 & loeduc == 1 & (runiform() <= .05)
replace good2 = 1 if good1 == 0 & loeduc == 0 & (runiform() <= .20)
replace good2 = 1 if good1 == 0 & treatment == 1 & (runiform() <= .30)
label val good2 health
// Transform some more and clean up
generate poor2 = (good2 == 0)
generate poor1 = (good1 == 0)
drop educ sex treatmentstr health1 number good1 good2
// Table 1-ish
table poor2, by(loeduc female treatment) contents(freq)
// 1) Regression adjustment
logit poor2 treatment female if loeduc == 0, or
logit poor2 treatment female if loeduc == 1, or
logit poor2 treatment##loeduc female##loeduc, or
// 2) Propensity score matching
teffects nnmatch (poor2 female) (treatment) if loeduc == 0
teffects nnmatch (poor2 female) (treatment) if loeduc == 1
// 3) Difference in difference
// Transform to long format
gen id = _n
reshape long poor, i(id) j(year)
logit poor treatment##c.year if loeduc == 0, or
logit poor treatment##c.year if loeduc == 1, or
logit poor loeduc##treatment##c.year, or
// 4) Fixed effects model
replace treatment = 0 if year == 1
xtset id year
xtreg poor treatment year if loeduc == 0, fe
xtreg poor treatment year if loeduc == 1, fe
xtreg poor treatment##loeduc year##loeduc, fe
Reference
Hu, Yannan, Frank J. van Lenthe, Rasmus Hoffmann, Karen van Hedel, and Johan P. Mackenbach. 2017. "Assessing the Impact of Natural Policy Experiments on Socioeconomic Inequalities in Health. How to Apply Commonly Used Quantitative Analytical Methods?" BMC Medical Research Methodology 17(1):68. doi: 10.1186/s12874-017-0317-5
Jul 29, 2017
OLS regression: overfitting and dichotomization
Babyak (2004) demonstrates a couple of aspects of OLS regression, making use of the following simulations:
// Figure 1
set seed 1
// Create file to store simulation results
tempname foo
postfile `foo' b using clt, replace
// Simulate analyses
forvalues i = 1/10000 {
drop _all
qui set obs 100
generate x = rnormal()
generate e = rnormal()
generate y = .4 * x + e
qui regress y x
local b = _b[x]
post `foo' (`b')
}
postclose `foo'
// Open results from simulations and plot
use clt, clear
histogram b, freq ytitle("Frequency of b value") xtitle("Values of b") name(figure1, replace)
// Figure 2
set seed 1
// Create file to store simulation results
tempname foo
postfile `foo' r2_50 r2_100 r2_150 r2_200 using overfitting, replace
// Simulate analyses
forvalues i = 1/10000 {
drop _all
set obs 10000
generate y = rnormal()
foreach i of numlist 1/15 {
generate x`i' = rnormal()
}
foreach j of numlist 50 100 150 200 {
preserve
sample `j', count
qui reg y x*
local r2_`j' = e(r2)
restore
}
post `foo' (`r2_50') (`r2_100') (`r2_150') (`r2_200')
}
postclose `foo'
// Open results from simulations
use overfitting, clear
// Plot
twoway (kdensity r2_200) ///
(kdensity r2_150) ///
(kdensity r2_100) ///
(kdensity r2_50) ///
, ytitle("Percent of samples") ///
xtitle("R-square value from regression model") ///
xlabel(0 (.1) .6) ///
ylabel(0 (2) 20) ///
legend(order(1 "ca. 13 cases/predictor ({it:N} = 200)" ///
2 "10 cases/predictor ({it:N} = 150)" ///
3 "ca. 7 cases/predictor ({it:N} = 100)" ///
4 "ca. 3 cases/predictor ({it:N} = 50)") ///
pos(2) ring(0)) name(figure2, replace)
// Figure 2
set seed 1
// Create file to store simulation results
tempname foo
postfile `foo' r2_50 r2_100 r2_150 r2_200 using overfitting, replace
// Simulate analyses
forvalues i = 1/10000 {
drop _all
set obs 10000
generate y = rnormal()
foreach i of numlist 1/15 {
generate x`i' = rnormal()
}
foreach j of numlist 50 100 150 200 {
preserve
sample `j', count
qui reg y x*
local r2_`j' = e(r2)
restore
}
post `foo' (`r2_50') (`r2_100') (`r2_150') (`r2_200')
}
postclose `foo'
// Open results from simulations
use overfitting, clear
// Plot
twoway (kdensity r2_200) ///
(kdensity r2_150) ///
(kdensity r2_100) ///
(kdensity r2_50) ///
, ytitle("Percent of samples") ///
xtitle("R-square value from regression model") ///
xlabel(0 (.1) .6) ///
ylabel(0 (2) 20) ///
legend(order(1 "ca. 13 cases/predictor ({it:N} = 200)" ///
2 "10 cases/predictor ({it:N} = 150)" ///
3 "ca. 7 cases/predictor ({it:N} = 100)" ///
4 "ca. 3 cases/predictor ({it:N} = 50)") ///
pos(2) ring(0)) name(figure2, replace)
// Figure 4 (actually Table 1)
clear
set seed 1
// Create file to store simulation results
tempname foo
postfile `foo' n correlation typei using dichotomization, replace
// Simulate analyses
forvalues i = 1/10000 {
drop _all
foreach j of numlist 50 100 200 {
foreach k of numlist 0 .3 .5 .7 {
qui drawnorm x1 x2, ///
n(`j') ///
corr(1, `k', 1) cstorage(lower) ///
clear
generate e = rnormal()
generate y = .5*x1 + 0*x2 + e
qui sum x1, detail
generate x1s = (x1 > r(p50))
qui sum x2, detail
generate x2s = (x2 > r(p50))
qui reg y x1s x2s
local typei = _b[x2s]/_se[x2s]
local sig = (abs(`typei') > 1.96)
*di _b[x2s] _skip(5) _se[x2s] _skip(5) `typei' _skip(5) `sig'
post `foo' (`j') (`k') (`sig')
}
}
}
postclose `foo'
// Open results from simulations
use dichotomization, clear
collapse typei, by(n correlation)
graph hbar typei, over(n, relabel(1 "{it:N} = 50" 2 "{it:N} = 100" 3 "{it:N} = 200")) ///
over(correlation, relabel(1 "{it:Corr(x{sub:1}, x{sub:2})} = 0" ///
2 "{it:Corr(x{sub:1}, x{sub:2})} = .3" ///
3 "{it:Corr(x{sub:1}, x{sub:2})} = .5" ///
4 "{it:Corr(x{sub:1}, x{sub:2})} = .7")) ///
ytitle("Type I error rate") yscale(alt) ylabel(, format(%6.2f)) ///
name(figure4, replace)
Reference
Babyak, Michael A. 2004. "What You See May Not Be What You Get. A Brief, Nontechnical Introduction to Overfitting in Regression-Type Models." Psychosomatic Medicine 66(3):411-421. doi: 10.1097/01.psy.0000127692.23278.a9
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
Jul 7, 2017
Analyzing censored data with the Tobit model
This replicates Table 2.2 in Breen (1996).
clear // Simulate data set obs 2000 generate ui = rnormal(0, 2) generate xi = rnormal() generate yi_star = 1 + 2*xi + ui drop ui // Fit OLS model regress yi_star xi eststo ols1 estadd scalar sigma = e(rmse) // Truncate variable generate yi = yi_star if yi_star > 0 replace yi = 0 if yi_star <= 0 recode yi (0 = .), gen(yi_h) // (A) OLS (using all observations on y // including y1 = 0) regress yi xi eststo ols2 estadd scalar sigma = e(rmse) // (B) OLS (yi > 0) regress yi xi if yi >0 eststo ols3 estadd scalar sigma = e(rmse) // (C) Heckman 2-step heckman yi_h xi, select(xi) twostep eststo heckman // (D) Tobit tobit yi xi, ll(0) eststo tobit estadd scalar sigma = _b[sigma:_cons] // Table 2.2esttab ols2 ols3 heckman tobit ols1, b(3) se(3) nostar stat(sigma) /// mtitles("(A) OLS incl. yi = 0" /// "(B) OLS yi > 0" /// "(C) Heckman 2-step" /// "(D) Tobit" /// "OLS yi_star") /// coeflabel(_cons "alpha" xi "beta") /// collabels() /// drop(mills:lambda sigma:_cons) /// order(_cons xi) /// unstack /// modelwidth(20) nonumber
Reference
Breen, Richard. 1996. Regression Models. Censored, Sample Selected, or Truncated Data. Sage. doi: 10.4135/9781412985611
Jul 4, 2017
Chapters 9 and 10 of Singer and Willett's (2003) book on longitudinal data analysis
// Figure 9.1use "C:\singer willett (2003)\teachers.dta", clear label define censor 0 "Not censored" 1 "Censored" label val censor censor histogram t, by(censor, legend(off) note("")) freq xlabel(1 (1) 12) ylabel(0 (100) 500) /// addlabels discrete xtitle("Years of teaching") name(figure91, replace) // Table 10.1generate event = !censor qui ltable t event, noadjust saving(lifetable, replace) preserve use lifetable, clear list t0 t1 start deaths lost hazard survival, sep(0) noobs erase lifetable.dta restore // Figure 10.1stset t, failure(event) sts generate h = h twoway scatter h t, msymbol(i) connect(l) ylabel(0 (.05) .15) xlabel(0 (1) 13) /// sort xtitle("Years in teaching") /// ytitle("Estimated hazard probability") /// name(figure101a, replace) nodraw ltable t event, noadjust notab graph noconf xlabel(0 (1) 13) ylabel(0 (.5) 1) /// yline(.5) xline(7.6) /// xtitle("Years in teaching") /// ytitle("Proportion surviving") /// name(figure101b, replace) nodraw graph combine figure101a figure101b, col(1) name(figure101, replace) ysize(8) // Figure 10.2use "C:\singer willett (2003)\relapse_days.dta", clear generate weeks = int(days / 7) + 1 generate event = !censor stset weeks, failure(event) sts generate h = h twoway (scatter h weeks, connect(l) msymbol(i)), xlabel(0 (1) 12) ylabel(, format(%6.2f)) /// ytitle("Estimated hazard probability") /// xtitle("Weeks after release") /// name(figure102a1, replace) nodraw ltable weeks event, notable graph noconf ylabel(0 (.25) 1, format(%6.2f)) /// xlabel(0(1)12) yline(0.5) /// ytitle("Estimated survival probability") /// xtitle("Weeks after release") /// name(figure102a2, replace) nodraw graph combine figure102a1 figure102a2, col(2) /// title("{bf:A} Time to cocaine relapse", pos(11)) name(figure102a, replace) nodraw use "C:\singer willett (2003)\firstsex.dta", clear generate event = !censor stset time, failure(event) sts generate h = h twoway (scatter h time, connect(l) msymbol(i) sort), xlabel(6 (1) 12) ylabel(, format(%6.2f)) /// ytitle("Estimated hazard probability") /// xtitle("Grade") /// name(figure102b1, replace) nodraw ltable time event, notable graph noconf ylabel(0 (.25) 1, format(%6.2f)) /// xlabel(6 (1) 12) yline(0.5) /// ytitle("Estimated survival probability") /// xtitle("Grade") /// name(figure102b2, replace) nodraw graph combine figure102b1 figure102b2, col(2) /// title("{bf:B} Age at first intercourse for males", pos(11)) name(figure102b, replace) nodraw use "C:\singer willett (2003)\suicide_orig.dta", clear generate event = !censor stset time, failure(event) sts generate h = h twoway (scatter h time, connect(l) msymbol(i) sort), xlabel(5 (2) 21) ylabel(, format(%6.2f)) /// ytitle("Estimated hazard probability") /// xtitle("Age") /// name(figure102c1, replace) nodraw ltable time event, notable graph noconf ylabel(0 (.25) 1, format(%6.2f)) /// xlabel(5 (2) 21) yline(0.5) /// ytitle("Estimated survival probability") /// xtitle("Age") /// name(figure102c2, replace) nodraw graph combine figure102c1 figure102c2, col(2) /// title("{bf:C} Age at first suicide ideation", pos(11)) name(figure102c, replace) nodraw use "C:\singer willett (2003)\congress_orig.dta", clear generate event = !censor stset time, failure(event) sts generate h = h twoway (scatter h time, connect(l) msymbol(i) sort), xlabel(0 (1) 8) ylabel(0 (.1) .3, format(%6.2f)) /// ytitle("Estimated hazard probability") /// xtitle("Terms in office") /// name(figure102d1, replace) nodraw ltable time event, notable graph noconf ylabel(0 (.25) 1, format(%6.2f)) /// xlabel(0 (1) 8) yline(0.5) /// ytitle("Estimated survival probability") /// xtitle("Terms in office") /// name(figure102d2, replace) nodraw graph combine figure102d1 figure102d2, col(2) /// title("{bf:D} Duration of congressional careers for females", pos(11)) name(figure102d, replace) nodraw graph combine figure102a figure102b figure102c figure102d, row(4) ysize(11) name(figure102, replace) // Table 10.2 use "C:\singer willett (2003)\teachers.dta", clear generate event = !censor qui ltable t event, noadjust saving(lifetable, replace) preserve use lifetable, clear list t0 hazard sehazard survival sesurvival, noobs sep(0) erase lifetable.dta restore // Figure 10.4 use "C:\singer willett (2003)\teachers.dta", clear list id t censor if inlist(id, 20, 126, 129), noobs sep(0) expand t bysort id: generate period = _n // Period identifier bysort id: generate event = !censor & _n == _N // Calculate outcome for // discrete-time analysis list id period event if inlist(id, 20, 126, 129), noobs sepby(id) // Table 10.3 table period, c(sum event n event mean event)
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
Chapter 5 of Singer and Willett's (2003) book on longitudinal data analysis
There is another take on the chapter here, but I like mine better.
// Table 5.1 use "C:\singer willett (2003)\reading_pp.dta", clear format age %6.2f list id wave agegrp age piat if inlist(id, 4, 27, 31, 33, 41, 49, 69, 77, 87), sep(0) noobs // Figure 5.1twoway (scatter piat age) /// (lfit piat age) /// (scatter piat agegrp) /// (lfit piat agegrp) /// if inlist(id, 4, 27, 31, 33, 41, 49, 69, 77, 87) /// , by(id, note("")) xtitle("{it:AGE} or {it:AGEGRP}") /// ylabel(0 (20) 80) xlabel(6 (1) 12, format(%6.0f)) ytitle("{it:PIAT}") /// legend(order(1 "Age" 2 "Linear fit age" /// 3 "Target age" 4 "Linear fit target age") col(2)) /// name(figure51, replace) ysize(8) // Table 5.2generate agegrp_65 = agegrp - 6.5 generate age_65 = age - 6.5 capture program drop randomslopetable program define randomslopetable eststo `1' qui estadd scalar dev = -2*e(ll) // Deviance qui matrix foo = e(N_g) // Number of individuals qui estadd scalar nc = foo[1,1] // Number of individuals qui estadd scalar v2 = exp(2*[lns1_1_1]_b[_cons]) // Slope variance qui estadd scalar v1 = exp(2*[lns1_1_2]_b[_cons]) // Intercept variance qui estadd scalar cov = tanh([atr1_1_1_2]_b[_cons]) * /// Slope-intercept covariance exp([lns1_1_1]_b[_cons]) * /// exp([lns1_1_2]_b[_cons]) qui estadd scalar v_e = exp(2*[lnsig_e]_b[_cons]) // Residual variance end eststo clear mixed piat agegrp_65 || id: agegrp_65, var cov(uns) mle eststo agegrp randomslopetable agegrp mixed piat age_65 || id: age_65, var cov(uns) mle eststo age randomslopetable age esttab /// , b(4) se(4) star(~ 0.10 * 0.05 ** 0.01 *** 0.001) /// // Re-define stars rename(agegrp_65 age_65) /// // Align coefficients coeflabels(_cons "Intercept" age_65 "Change") /// // Label coefficients order(_cons age_65) /// // Order coefficients stats(v_e v1 v2 dev aic bic nc N, /// Add variance components to table fmt(2 2 2 1 1 1 0 0) /// labels("Var(Residual)" /// "Var(Initial)" /// "Var(Change)" /// "Deviance" /// "AIC" /// "BIC" /// "No. individuals" /// "No. measurements")) /// nonumbers nobaselevels noomitted varwidth(25) /// mtitles("AGEGRP - 6.5" "AGE - 6.5") /// keep(piat:) /// // Drop variance components in weird shapes eqlabels("", none) // Removes equation label // Table 5.3 use "C:\singer willett (2003)\wages_pp.dta", clear list id exper lnw black hgc uerate if inlist(id, 206, 332, 1028), sep(0) noobs // Table 5.4eststo clear mixed lnw exper || id: exper, var mle cov(uns) eststo modela randomslopetable modela mixed lnw c.exper##c.hgc_9 c.exper##black || id: exper, var mle cov(uns) eststo modelb randomslopetable modelb mixed lnw exper hgc_9 c.exper#black || id: exper, var mle cov(uns) eststo modelc randomslopetable modelc esttab /// , b(4) se(4) star(~ 0.10 * 0.05 ** 0.01 *** 0.001) /// // Re-define stars coeflabels(_cons "Intercept" hgc_9 "(HGC - 9)" /// 1.black "BLACK" exper "Change" /// c.exper#c.hgc_9 "Change x (HGC - 9)" /// 1.black#c.exper "Change x BLACK") /// // Label coefficients order(_cons hgc_9 1.black /// exper c.exper#c.hgc_9 1.black#c.exper) /// // Order coefficients stats(v_e v1 v2 dev aic bic nc N, /// Add variance components to table fmt(4 4 4 1 1 1 0 0) /// labels("Var(Residual)" /// "Var(Initial)" /// "Var(Change)" /// "Deviance" /// "AIC" /// "BIC" /// "No. individuals" /// "No. measurements")) /// nonumbers nobaselevels noomitted varwidth(25) /// mtitles("Model A" "Model B" "Model C") /// keep(lnw:) /// // Drop variance components in weird shapes eqlabels("", none) // Removes equation label // Figure 5.2estimates restore modelc margins, at(exper = (0 (1) 10) black = (0 1) hgc_9 = (0 3)) marginsplot, ytitle("Predicted log hourly wage") title("") /// recastci(rarea) ciopt(color(gs14)) /// xtitle("{it:EXPER}") /// addplot(scatteri 2.35 10 "White", msymbol(none) mlabpos(0) /// || scatteri 2.23 10 "White", msymbol(none) mlabpos(0) /// || scatteri 2.185 10 "Black", msymbol(none) mlabpos(0) /// || scatteri 2.07 10 "Black", msymbol(none) mlabpos(0) /// || scatteri 1.75 0 "9th grade", msymbol(none) mlabpos(0) /// || scatteri 1.88 0 "12th grade", msymbol(none) mlabpos(0)) /// legend(off) ylabel(, format(%6.1f)) /// name(figure52, replace) // Table 5.5use "C:\singer willett (2003)\wages_small_pp.dta", clear eststo clear mixed lnw hgc_9 exper c.exper#black || id: exper, var mle cov(uns) eststo modela randomslopetable modela mixed lnw hgc_9 exper c.exper#black || id: , var mle cov(uns) eststo modelc qui estadd scalar dev = -2*e(ll) // Deviance qui matrix foo = e(N_g) // Number of individuals qui estadd scalar nc = foo[1,1] // Number of individuals qui estadd scalar v1 = exp(2*[lns1_1_1]_b[_cons]) // Intercept variance qui estadd scalar v_e = exp(2*[lnsig_e]_b[_cons]) // Residual variance esttab /// , b(4) se(4) star(~ 0.10 * 0.05 ** 0.01 *** 0.001) /// // Re-define stars rename(agegrp_65 age_65) /// // Align coefficients coeflabels(_cons "Intercept" hgc_9 "(HGC - 9)" /// 1.black "BLACK" exper "Change" /// 1.black#c.exper "Change x BLACK") /// // Label coefficients order(_cons hgc_9 1.black /// exper 1.black#c.exper) /// // Order coefficients stats(v_e v1 v2 dev aic bic nc N, /// Add variance components to table fmt(4 4 4 1 1 1 0 0) /// labels("Var(Residual)" /// "Var(Initial)" /// "Var(Change)" /// "Deviance" /// "AIC" /// "BIC" /// "No. individuals" /// "No. measurements")) /// nonumbers nobaselevels noomitted varwidth(25) /// mtitles("Model A" "Model C") /// keep(lnw:) /// // Drop variance components in weird shapes eqlabels("", none) // Removes equation label // Table 5.6 use "C:\singer willett (2003)\unemployment_pp.dta", clear list id months cesd unemp if inlist(id, 7589, 55697, 67641, 65441, 53782), sepby(id) noobs // Table 5.7eststo clear mixed cesd months || id: months, var mle cov(uns) eststo modela randomslopetable modela mixed cesd months i.unemp || id: months, var mle cov(uns) eststo modelb randomslopetable modelb mixed cesd c.months##i.unemp || id: months, var mle cov(uns) eststo modelc randomslopetable modelc version 10 // This one is a challenge to fit generate unempXmonths = unemp * months xtmixed cesd unemp unempXmonths || id: unemp unempXmonths, var mle cov(uns) eststo modeld qui estadd scalar dev = -2*e(ll) // Deviance qui matrix foo = e(N_g) // Number of individuals qui estadd scalar nc = foo[1,1] // Number of individuals qui estadd scalar v1 = exp(2*[lns1_1_1]_b[_cons]) // Intercept variance qui estadd scalar v4 = exp(2*[lns1_1_2]_b[_cons]) // Unemp variance qui estadd scalar v3 = exp(2*[lns1_1_3]_b[_cons]) // Unemp x months variance qui estadd scalar v_e = exp(2*[lnsig_e]_b[_cons]) // Residual variance version 14 esttab /// , b(4) se(4) star(~ 0.10 * 0.05 ** 0.01 *** 0.001) /// // Re-define stars rename(unemp 1.unemp unempXmonths 1.unemp#c.months) /// coeflabels(_cons "Intercept" months "Change" /// 1.unemp "UNEMP" /// 1.unemp#c.months "Change x UNEMP") /// // Label coefficients order(_cons months /// 1.unemp 1.unemp#c.months) /// // Order coefficients stats(v_e v1 v2 v3 v4 dev aic bic nc N, /// Add variance components to table fmt(4 4 4 4 4 1 1 1 0 0) /// labels("Var(Residual)" /// "Var(Initial)" /// "Var(Change)" /// "Var(UNEMP)" /// "Var(UNEMP x TIME)" /// "Deviance" /// "AIC" /// "BIC" /// "No. individuals" /// "No. measurements")) /// nonumbers nobaselevels noomitted varwidth(25) /// mtitles("Model A" "Model B" "Model C" "Model D") /// keep(cesd:) /// // Drop variance components in weird shapes eqlabels("", none) // Removes equation label // Figure 5.3estimates restore modelb margins, at(months = (0 15) unemp = (0 1)) marginsplot, ylabel(5 (5) 20) /// recastci(rarea) ciopts(color(gs14)) /// ytitle("{it:Predicted CES-D}") /// xtitle("Months since job loss") /// addplot(scatteri 11 15 "Employed", msymbol(none) mlabpos(9) /// || scatteri 16 15 "Unemployed", msymbol(none) mlabpos(9) /// xlabel(0 (2) 14)) /// legend(off) title("") /// name(figure53, replace) ** All other plots of Figure 5.3 are only variants of this one ** // Figure 5.4 estimates restore modelb margins, at(months = (0 (1) 15) unemp = (0 1)) marginsplot, ylabel(5 (5) 20) /// recastci(rarea) ciopts(color(gs14)) /// ytitle("Predicted {it:CES-D}") /// xtitle("Months since job loss") /// addplot(scatteri 11 15 "Employed", msymbol(none) mlabpos(9) /// || scatteri 16 15 "Unemployed", msymbol(none) mlabpos(9) /// xlabel(0 (2) 14)) /// legend(off) title("Model B") subtitle("Main effects of" /// "{it:UNEMP} and {it:TIME}") /// name(figure54A, replace) estimates restore modelc margins, at(months = (0 (1) 15) unemp = (0 1)) marginsplot, ylabel(5 (5) 20) /// recastci(rarea) ciopts(color(gs14)) /// ytitle("Predicted {it:CES-D}") /// xtitle("Months since job loss") /// addplot(scatteri 11 15 "Employed", msymbol(none) mlabpos(9) /// || scatteri 16 15 "Unemployed", msymbol(none) mlabpos(9) /// xlabel(0 (2) 14)) /// legend(off) title("Model C") subtitle("Interaction between" /// "{it:UNEMP} and {it:TIME}") /// name(figure54B, replace) estimates restore modeld // not sure how to do this with marginsplot predict pd twoway (line pd months if unemp == 0, c(L)) /// (line pd months if unemp == 1, c(L)) /// (scatteri 11.5 15 "Employed", msymbol(none) mlabpos(9)) /// (scatteri 15 15 "Unemployed", msymbol(none) mlabpos(9)) /// , ylabel(5 (5) 20) legend(off) /// xlabel(0 (2) 14) /// ytitle("Predicted {it:CES-D}") /// xtitle("Months since job loss") /// title("Model D") subtitle("Constraining the effect of {it:TIME}" /// "among the re-employed") /// name(figure54C, replace) graph combine figure54A figure54B figure54C, row(1) xsize(11) // Table 5.8 use "C:\singer willett (2003)\wages_pp.dta", clear eststo clear mixed lnw hgc_9 ue_7 exper c.exper#black || id: exper, mle cov(un) var eststo modela randomslopetable modela mixed lnw hgc_9 ue_mean ue_person_centered exper c.exper#black || id: exper, mle cov(un) var eststo modelb randomslopetable modelb mixed lnw hgc_9 ue1 ue_centert1 exper c.exper#black || id: exper, mle cov(un) var eststo modelc randomslopetable modelc esttab /// , b(4) se(4) star(~ 0.10 * 0.05 ** 0.01 *** 0.001) /// // Re-define stars rename(ue_mean ue_7 ue1 ue_7 ue_centert1 ue_person_centered) /// coeflabels(_cons "Intercept" hgc_9 "(HGC - 9)" /// ue_7 "UERATE" ue_person_centered "Deviation UERATE" /// exper "Change" 1.exper#black "Change x BLACK") /// // Label coefficients order(_cons hgc_9 ue_7 ue_person_centered /// exper c.exper#black) /// // Order coefficients stats(v_e v1 v2 dev aic bic nc N, /// Add variance components to table fmt(4 4 4 1 1 1 0 0) /// labels("Var(Residual)" /// "Var(Initial)" /// "Var(Change)" /// "Deviance" /// "AIC" /// "BIC" /// "No. individuals" /// "No. measurements")) /// nonumbers nobaselevels noomitted varwidth(25) /// mtitles("Model A" "Model B" "Model C") /// keep(lnw:) /// // Drop variance components in weird shapes eqlabels("", none) // Removes equation label // Table 5.9 use "C:\singer willett (2003)\medication_pp.dta", clear list wave day timeofday time time333 time667 in 1/11, noobs sep(0) // Table 5.10 eststo clear mixed pos i.treat##c.time || id: time, mle cov(uns) var eststo modela randomslopetable modela mixed pos i.treat##c.time333 || id: time333, mle cov(uns) var eststo modelb randomslopetable modelb mixed pos i.treat##c.time667 || id: time667, mle cov(uns) var eststo modelc randomslopetable modelc esttab /// , b(2) se(2) star(~ 0.10 * 0.05 ** 0.01 *** 0.001) /// // Re-define stars rename(time333 time time667 time /// 1.treat#c.time333 1.treat#c.time /// 1.treat#c.time667 1.treat#c.time) /// coeflabels(_cons "Intercept" time "Change" /// 1.treat "TREAT" 1.treat#c.time "Change x TREAT") /// // Label coefficients order(_cons 1.treat time 1.treat#c.time) /// // Order coefficients stats(v_e v1 v2 cov dev aic bic nc N, /// Add variance components to table fmt(2 2 2 2 1 1 1 0 0) /// labels("Var(Residual)" /// "Var(Initial)" /// "Var(Change)" /// "Cov(Init., Change)" /// "Deviance" /// "AIC" /// "BIC" /// "No. individuals" /// "No. measurements")) /// nonumbers nobaselevels noomitted varwidth(25) /// mtitles("Model A" "Model B" "Model C") /// keep(pos:) /// // Drop variance components in weird shapes eqlabels("", none) // Removes equation label // Figure 5.5estimates restore modela margins, at(time = (0 (1) 7) treat = (1 0)) marginsplot, recastci(rarea) ciopts(color(gs14)) legend(off) /// title("") xtitle("Days") ytitle("Predicted {it:POS}") /// addplot(scatteri 187 7 "Treatment", msymbol(none) mlabpos(11) /// || scatteri 153 7 "Control", msymbol(none) mlabpos(11) /// xlabel(0 (1) 7)) /// name(figure55, 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
Subscribe to:
Posts (Atom)