Reference Mummolo, Jonathan, and Erik Peterson. 2018. "Improving the Interpretation of Fixed Effects Regression Results." Political Science Research and Methods 6(4):829-835. doi: 10.1017/psrm.2017.44
set scheme lean1 webuse nlswork, clear xtset idcode year // Model of interest: xtreg ln_w ttl_exp age c.age#c.age tenure c.tenure#c.tenure not_smsa south i.year, fe local b : display %4.2f _b[ttl_exp] // Isolate relevant variation in the treatment: reghdfe ttl_exp, absorb(idcode year) residuals(relevantvariation) // Identify a plausible counterfactual shift in X given the data: center ttl_exp sum relevantvariation local sd1 : display %4.2f r(sd) sum ttl_exp local sd2 : display %4.2f r(sd) twoway (kdensity c_ttl_exp) /// (kdensity relevantvariation), /// ytitle("Density") xtitle(Treatment (centered around 0)) /// title("{bf:A} Treatment distribution", span bexpand justification(left)) /// legend(order(1 "Original distribution, {it:SD} = `sd2'" /// 2 "Distribution after FE, {it:SD} = `sd1'" ) /// pos(2) ring(0)) /// name(figure1A, replace) // Calculate range within individuals bysort idcode: egen max = max(ttl_exp) bysort idcode: egen min = min(ttl_exp) generate range = max - min drop max min qui sum range, detail local mean : display %4.2f r(mean) local p95 : display %4.2f r(p95) twoway (histogram range, freq) /// (scatteri 2000 `mean' (3) "Mean = `mean'", msymbol(i)) /// (scatteri 2000 `p95' (3) "95th percentile = `p95'" , msymbol(i)), /// title("{bf:B} Within-individual range of treatment", span bexpand justification(left)) /// ylabel(, format(%7.0gc)) ytitle(Frequency) /// xline(`mean') xline(`p95') /// xtitle(Within-individual range) legend(off) /// name(figure1B, replace) graph combine figure1A figure1B, col(2) // Calculate % of cases where treatment does not vary egen pickone = tag(idcode) fre range if pickone // In 12% of cases the treatment does not vary // Multiply the estimated coefficient of interest by the revised standard deviation display "Estimated coefficient of interest = `b'" display "Within-standard deviation = `sd1'" display "Product = " _skip(1) `sd1' * `b'
Showing posts with label twoway kdensity. Show all posts
Showing posts with label twoway kdensity. Show all posts
May 8, 2024
Mummolo & Peterson (2018): Fixed-effects regression
May 16, 2018
Random graphs (132): Distributions
// Open ESS 2006 data
use ESS3e03_6.dta, clear
// Prepare variables
kountry cntry, from(iso2c)
rename NAMES_STD country
recode agea (999 = .)
recode ygcdbyr (6666/9999 = .)
recode inwyye (9999 = .)
recode gndr (9 = .)
// Generate relevant variables
generate yearborn = inwyye - agea
gen grandparent = ygcdbyr - yearborn
//Plot figure
twoway kdensity grandparent if gndr == 1, by(country, title(Men , justification(left) span) note("")) xline(50) ylabel(0 .05 .10, format(%6.2f)) ytitle("Density") xtitle("") /*xtitle("Age at grandparenthood")*/ name(men, replace) nodraw
twoway kdensity grandparent if gndr == 2, by(country, title(Women, justification(left) span) note("")) xline(50) ylabel(0 .05 .10, format(%6.2f)) ytitle("Density") xtitle("Age at grandparenthood") name(women, replace) nodraw
graph combine men women, col(1) ysize(8) note("{it:Note:} Vertical line indicates age 50", size(vsmall))
Mar 9, 2018
Random graphs (127): Density plots
clear
input str25 inout days
"Decision" 15
"Decision" 15
"Decision" 25
"Decision" 1
"Decision" 135
"Decision" 12
"Decision" 184
"Decision" 40
"Decision" 12
"Decision" 11
"Decision" 44
"Decision" 35
"Decision" 87
"Decision" 101
"Decision" 196
"Decision" 51
"Decision" 120
"Decision" 7
"Decision" 84
"Decision" 103
"Decision" 2
"Decision" 175
"Decision" 29
"Decision" 384
"Decision" 56
"Decision" 49
"Decision" 28
"Decision" 40
"Decision" 140
"Decision" 3
"Decision" 50
"Decision" 5
"Decision" 85
"Decision" 43
"Decision" 152
"Decision" 19
"Decision" 11
"Decision" 103
"Decision" 27
"Decision" 133
"Decision" 64
"Decision" 68
"Decision" 113
"Decision" 159
"Decision" 19
"Decision" 82
"Decision" 44
"Decision" 2
"Decision" 75
"Decision" 1
"Decision" 90
"Decision" 10
"Decision" 110
"Reviewed" 21
"Reviewed" 1
"Declined" 0
"Reviewed" 30
"Reviewed" 0
"Declined" 0
"Reviewed" 0
"Declined" 0
"Declined" 1
"Reviewed" 6
"Declined" 0
"Declined" 0
"Declined" 1
"Declined" 1
"Reviewed" 47
"Declined" 0
"Declined" 0
"Declined" 1
"Reviewed" 1
"Declined" 0
"Reviewed" 2
"Reviewed" 9
"Declined" 2
"Reviewed" 1
"Reviewed" 7
"Reviewed" 28
"Reviewed" 2
"Declined" 0
"Reviewed" 16
"Reviewed" 89
"Reviewed" 4
"Declined" 1
"Reviewed" 3
"Declined" 4
"Reviewed" 38
"Reviewed" 1
"Reviewed" 41
"Reviewed" 28
"Reviewed" 51
"Declined" 0
"Reviewed" 10
"Reviewed" 59
"Declined" 1
"Reviewed" 13
"Declined" 0
"Reviewed" 6
"Reviewed" 12
"Reviewed" 45
"Reviewed" 25
"Reviewed" 4
"Reviewed" 8
"Reviewed" 0
"Reviewed" 3
"Reviewed" 7
"Reviewed" 64
"Reviewed" 49
"Reviewed" 4
"Reviewed" 57
"Reviewed" 1
"Reviewed" 1
"Reviewed" 0
"Reviewed" 1
"Reviewed" 42
"Reviewed" 10
"Declined" 0
"Reviewed" 49
"Reviewed" 4
"Reviewed" 3
"Declined" 0
"Reviewed" 16
"Reviewed" 46
"Reviewed" 37
"Reviewed" 2
"Reviewed" 14
"Reviewed" 7
"Reviewed" 21
"Reviewed" 3
"Reviewed" 34
"Reviewed" 14
"Reviewed" 4
"Reviewed" 29
"Reviewed" 20
end
// Using label rather than order for the legend labels allows using a line break in the legend
twoway (kdensity days if inout == "Decision", lcolor(red)) ///
(kdensity days if inout == "Reviewed" | inout == "Declined", lcolor(green) lpattern(solid)) ///
(kdensity days if inout == "Reviewed", lcolor(midgreen) lpattern(solid)), ///
legend(label(1 "Waiting times to receive journal decisions") ///
label(2 "Waiting times to respond to review request") ///
label(3 "Waiting times to respond to review request" "(excluding when I declined to review)") ///
pos(2) ring(0)) ///
xlabel(0 25 50 75 100 200 300 400) xtitle(Days) ///
ytitle("Kernel density estimate") ///
title("How long journals wait for my reviews of their manuscripts" "{it:vs.} how long I wait for journals' decisions on my manuscripts", span)
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
Subscribe to:
Posts (Atom)