clear
// Generate data
set seed 1
set obs 50
gen hours = rnormal(3, 1) // Number of hours studied
gen e = rnormal(1,1)
gen questions = 2 + 2*hours + 1*e // Questions answered correctly
qui sum questions, detail
generate pass = (questions >= r(p75)) // Passing the exam
// 1) Histogram of outcome variable
twoway (histogram pass, discrete percent), ///
xlabel(0 "[0] Failed" 1 "[1] Passed") xtitle("") ///
ytitle("Percent of students") xsize(4) ysize(4) name(figure5, replace)
// 2) Scatterplot
twoway (scatter pass hours), ///
xlabel(0 (1) 5) xtitle("Hours studied for exam") ///
ytitle("Exam success") ///
ylabel(0 "[0] Failed" 1 "[1] Passed") legend(off) ///
xsize(4) ysize(4) name(figure6, replace)
// 3) Scatterplot with regression line
regress pass hours
local intercept = round(_b[_cons], .01)
local x = round(_b[hours], .01)
twoway (scatter pass hours) ///
(lfit pass hours, lpattern(solid) range(1 5)), ///
xlabel(0 (1) 5) xtitle("Hours studied for exam") ///
text(.8 2 "y = `intercept' + `x' x + e", size(large)) ///
ytitle("Exam success") ///
ylabel(0 "[0] Failed" 1 "[1] Passed") legend(off) ///
xsize(4) ysize(4) name(figure7, replace)
// 4) Logit curve
logit pass hours
predict yhat
twoway (scatter pass hours) ///
(line yhat hours, lpattern(solid) sort), ///
xlabel(0 (1) 5) xtitle("Hours studied for exam") ///
ytitle("Exam success") ///
ylabel(0 "[0] Failed" 1 "[1] Passed") legend(off) ///
xsize(4) ysize(4) name(figure8, replace)
graph combine figure5 figure6 figure7 figure8, ///
col(2) xsize(8) ysize(8) altshrink name(figures58, replace)
Showing posts with label twoway lfit. Show all posts
Showing posts with label twoway lfit. Show all posts
Jul 24, 2018
Random graphs (137): Logistic regression
Labels:
logit,
predict,
Random graphs,
Simulation,
twoway histogram,
twoway lfit,
twoway line,
twoway scatter
Random graphs (135): Scatterplot with OLS regression
clear
// Generate data
set seed 1
set obs 50
gen hours = rnormal(3, 1) // Number of hours studied
gen e = rnormal(1,1)
gen questions = 2 + 2*hours + 1*e // Questions answered correctly
// 1) Basic scatterplot
twoway (scatter questions hours), ///
xlabel(0 (1) 5) xtitle("Hours studied for exam") ///
ytitle("Number of questions answered correctly") ///
name(figure1, replace) ylabel(0 (5) 15) xsize(4) ysize(4)
// 2) Scatterplot with regression line
twoway (scatter questions hours) ///
(lfit questions hours), ///
xlabel(0 (1) 5) xtitle("Hours studied for exam") ///
ytitle("Number of questions answered correctly") ///
ylabel(0 (5) 15) legend(off) name(figure2, replace) xsize(4) ysize(4)
// 3) Scatterplot with regression line and equation
regress questions hours
local intercept = round(_b[_cons], .1)
local x = round(_b[hours], .1)
twoway (scatter questions hours) ///
(lfit questions hours), ///
xlabel(0 (1) 5) xtitle("Hours studied for exam") ///
text(14 2 "y = `intercept' + `x' + e", size(large)) ///
ytitle("Number of questions answered correctly") ///
ylabel(0 (5) 15) legend(off) name(figure3, replace) xsize(4) ysize(4)
// 4) Scatterplot with regression line and equation and labels for components
twoway (scatter questions hours) ///
(lfit questions hours, lpattern(solid) range(0 5)) ///
(function y = 1.6 + 2.5, range(1 2)) ///
(function y = 2, range(4.1 6.6) horizontal lpattern(solid) lcolor(red)) ///
(scatteri 1.6 0 (3) "Intercept", msymbol(o) mlabcolor(red)) ///
(scatteri 5 2 "Slope", msymbol(i) mlabcolor(red)), ///
xlabel(0 (1) 5) xtitle("Hours studied for exam") ///
text(14 2 "y = `intercept' + `x' x + e", size(large)) ///
ytitle("Number of questions answered correctly") ///
ylabel(0 1.6 5 10 15) legend(off) name(figure4, replace) xsize(4) ysize(4)
graph combine figure1 figure2 figure3 figure4, col(2) xsize(8) ysize(8) altshrink name(figures14, replace)
Jun 19, 2017
Random graphs (99): Scatterplots
regress treatments index
local r2 = round(e(r2), .01) // Round coeff.
twoway (scatter treatments index) ///
(lfit treatments index) ///
(scatter treatments index if inlist(cntry, "IT", "AT", "UA", "BE", "DK", "RU"), msymbol(o)) ///
, xtitle(ART comprehensive availability) ytitle("ART treatments" "per million women 15-44 y.") ///
title("{bf:B}", justification(left) bexpand span) ///
xlabel(0 (1) 9) ///
text(1500 6.6 "Denmark") ///
text(1500 8.6 "Belgium") ///
text(470 0.4 "Italy") ///
text(400 1.6 "Austria") ///
text(120 7.5 "Ukraine") ///
text(210 8.0 "Kazakhstan", placement(west)) ///
text(150 8.0 "Russia", placement(east)) ///
legend(order(2 "Linear fit R{char 178} = `r2'") ring(0) pos(7)) name(policy, replace)
Labels:
local,
Random graphs,
round(),
twoway lfit,
twoway scatter,
twoway scatteri
Dec 17, 2014
Random graphs (39): Plotting on a log axis
Rather than looking at a variable in absolute terms, namely GDP per capita in the upper panel of the Figure, it is also sometimes helpful to look at it in terms of percentage increases, as in the lower panel. Equal distances on the x-axis refer to equal percentage increases in GDP per capita. Four equidistant points on the x-axis are labeled, each indicating a fourfold increase in GDP.
use wvs2005_v20090901a.dta, clear
// Country variable
// 1) Turn into string
decode v2, gen(ctry)
// 2) Abbreviate
kountry ctry, from(other) stuck marker
ren _ISO3N_ cntry
kountry cntry, from(iso3n) to(iso2c)
ren _ISO2C_ country
// Generate outcome: % in good health
generate goodhealth = 100 if v11 <= 2
replace goodhealth = 0 if v11 > 2
replace goodhealth = . if v11 == .
// Collapse data set
collapse (mean) goodhealth [pw = v259], by(country)
// Generate year variable
gen year = 2006
// Preserve collapsed data set
preserve
// Get GDP from World Bank data base
wbopendata, language(en - English) country() topics() indicator(NY.GDP.PCAP.PP.CD) clear long
// Fix obtained data set
ren ny_gdp_pcap_pp_cd gdp
ren iso2code country
keep if year == 2006
keep gdp country
drop if country == ""
// Save obtained GDP data
tempfile gdp
save `gdp'
// Restore
restore
// Merge GDP with collapsed data set
merge m:1 country using `gdp'
keep if _merge == 3
drop _merge
// Create labels with thousand separator
label define gdp 20000 "20,000" ///
40000 "40,000" ///
60000 "60,000"
label val gdp gdp
// Plot on unlogged axis
twoway (scatter goodhealth gdp, mlabel(country) mlabpos(0) msymbol(none)) ///
(lfit goodhealth gdp) ///
, legend(off) xtitle("GDP per capita, 2006, PPP in current international $") ///
ytitle("% in good health") ///
xlabel(0(20000)60000, valuelabels) ///
name(unlogged, replace)
// Generate logged variable
gen loggdp = log(gdp) * 1000 // Multiply by 1,000 because only integers can be labeled
// Generate numbers for labeling
// Round them to three decimal digits, then multiply by 1,000 to get integers
local log1 = round(log(1000), .001) * 1000
local log2 = round(log(1000 * 4), .001) * 1000
local log3 = round(log(1000 * 4 * 4), .001) * 1000
local log4 = round(log(1000 * 4 * 4 * 4), .001) *1000
*di `log1' _skip(2) `log2' _skip(2) `log3' _skip(2) `log4'
// Create labels for numbers to be labeled
label define loggdp `log1' "1,000" ///
`log2' "4,000" ///
`log3' "16,000" ///
`log4' "64,000"
label value loggdp loggdp
// Plot on log axis
twoway (scatter goodhealth loggdp, mlabel(country) mlabpos(0) msymbol(none)) ///
(lfit goodhealth loggdp) ///
, legend(off) xtitle("GDP per capita, 2006, PPP in current international $") ///
ytitle("% in good health") ///
xlabel(`log1' `log2' `log3' `log4', valuelabels) ///
name(logged, replace)
// Combine plots
graph combine unlogged logged, col(1) ysize(8)
Labels:
_skip(#),
collapse,
decode,
graph combine,
kountry,
log(),
mlabel,
preserve,
Random graphs,
round(),
tempfile,
twoway lfit,
twoway scatter,
wbopendata,
World Values Survey,
ysize
Oct 24, 2014
Random graphs (34): Country scatterplots
use ART.dta, replace
// Fix missing values
recode avgemtrans (-1 -3 = .)
recode emhum (-1 = .)
// Calculate r-squared
qui regress avgemtrans emhum
local r2 = round(e(r2), .01) // Round coeff.
// Create first panel
twoway (scatter avgemtrans emhum, mlab(countrycode) mlabpos(0) msymbol(i)) ///
(lfit avgemtrans emhum) ///
, xlabel(, format(%6.1f)) ylabel(, format(%6.1f)) ///
ytitle("Avg. no. of fresh non-donor embryos" "per fresh embryo transfer cycle") ///
xtitle("Avg. level disagreement that" "embryo is a human being") ///
legend(order(2 "Linear fit R{char 178} = `r2'") ring(0) pos(8)) ///
name(figure2a, replace)
// Fix missing values
recode afford clonehelpinf (-1 = .)
// Listwise deletion to control y-axis range
preserve
keep if !missing(afford, clonehelpinf)
// Calculate r-squared
regress afford clonehelpinf
local r2 = round(e(r2), .01) // Round coeff.
// Create second panel
twoway (scatter afford clonehelpinf, mlab(countrycode) mlabpos(0) msymbol(i)) ///
(lfit afford clonehelpinf) ///
, xlabel(, format(%6.1f)) ylabel(0 (5) 25, format(%6.1f)) ///
ytitle("Affordability (net cost of a fresh ART cycle" "as % of annual disposable income)") ///
xtitle("Avg. level disagreement with" "cloning to help infertile couples") ///
legend(order(2 "Linear fit R{char 178} = `r2'") ring(0) pos(8)) ///
name(figure2b, replace)
restore
// Create final Figure
graph combine figure2a figure2b, row(1)
Jul 1, 2014
Random graphs (23): Between and within regression
clear
input j i xij x_j yij y_j
// Data from Snijders and Bosker (1999, Table 3.2):
1 1 1 2 5 6
1 2 3 2 7 6
2 1 2 3 4 5
2 2 4 3 6 5
3 1 3 4 3 4
3 2 5 4 5 4
4 1 4 5 2 3
4 2 6 5 4 3
5 1 5 6 1 2
5 2 7 6 3 2
end
// Plot as Figure 3.4 of Snijders ans Bosker (1999):
twoway (scatter yij xij) ///
(lfit yij xij) ///
(lfit y_j x_j) ///
(lfit yij xij if j == 1, lpattern(dash_dot)) ///
(lfit yij xij if j == 2, lpattern(dash_dot)) ///
(lfit yij xij if j == 3, lpattern(dash_dot)) ///
(lfit yij xij if j == 4, lpattern(dash_dot)) ///
(lfit yij xij if j == 5, lpattern(dash_dot)) ///
, legend(label(2 "Total regression") ///
label(3 "Between regression") ///
label(4 "Within regression") ///
order(2 3 4) ///
pos(1) ring(0)) ///
xtitle(X) ytitle(Y) ///
title("Within, between, and total relations") ///
xlabel(none) ylabel(none) ///
name(wbt_legend, replace)
// Lines directly labeled
twoway (scatter yij xij) ///
(lfit yij xij) ///
(lfit y_j x_j) ///
(lfit yij xij if j == 1, lpattern(dash_dot)) ///
(lfit yij xij if j == 2, lpattern(dash_dot)) ///
(lfit yij xij if j == 3, lpattern(dash_dot)) ///
(lfit yij xij if j == 4, lpattern(dash_dot)) ///
(lfit yij xij if j == 5, lpattern(dash_dot)) ///
, legend(off) ///
text(1.8 6.5 "Between") ///
text(3.0 7.5 "Total") ///
text(7.0 3.5 "Within") ///
text(5.0 5.5 "Within") ///
xtitle(X) ytitle(Y) ///
title("Within, between, and total relations") ///
xlabel(none) ylabel(none) ///
name(wbt_nolegend, replace)
graph combine wbt_legend wbt_nolegend, col(1) ysize(8) xcommon ycommon
// ysize(4) and xsize(5.5) are defaults
Reference
Snijders, Tom, and Roel Boskers. 1999. Multilevel Analysis. An Introduction to Basic and Advanced Multilevel Modeling. Sage.Jul 9, 2013
Morgan and Winship's (2007) example for bias due to conditioning on a collider in Stata
Morgan and Winship (2007: p. 66) illustrate Pearl's (2009) concern about conditioning on a collider variable using a simple example.
The general problem of conditioning on a collider is as follows. Consider three variables A, B, and C, with both A and B being causes of C: A → C ← B. (Formally, any variable C that has two arrows pointing to it along a given path is a collider.) Unlike a confounder (an uncontrolled common cause of A and B), a collider does not induce a zero-order correlation between A and B. However, when handled inappropriately, a collider can induce a conditional correlation between A and B.
Morgan and Winship's example shows just that: A college admits applicants based on their SAT scores and ratings of their motivation based on an interview. Those in the top 15 per cent of the sum of SAT and motivation ratings are being admitted. SAT scores and motivation ratings are largely uncorrelated.
What the Figures show is that the very small correlation between motivation and SAT score for the overall group turns out to be much larger when conditioning for admission status.
This also shows in an OLS regression:
Cole, Stephen R., Robert W. Platt, Enrique F. Schisterman, Haitao Chu, Daniel Westreich, David Richardson, and Charles Poole. 2010. "Illustrating Bias Due to Conditioning on a Collider." International Journal of Epidemiology 39(2):417-420. doi: 10.1093/ije/dyp334
Morgan, Stephen L., and Christopher Winship. 2007. Counterfactuals and Causal Inference. Methods and Principles for Social Research. Cambridge University Press.
Pearl, Judea. 2009. Causality. Models, Reasoning, and Inference, 2nd ed. Cambridge University Press.
The general problem of conditioning on a collider is as follows. Consider three variables A, B, and C, with both A and B being causes of C: A → C ← B. (Formally, any variable C that has two arrows pointing to it along a given path is a collider.) Unlike a confounder (an uncontrolled common cause of A and B), a collider does not induce a zero-order correlation between A and B. However, when handled inappropriately, a collider can induce a conditional correlation between A and B.
Morgan and Winship's example shows just that: A college admits applicants based on their SAT scores and ratings of their motivation based on an interview. Those in the top 15 per cent of the sum of SAT and motivation ratings are being admitted. SAT scores and motivation ratings are largely uncorrelated.
// Generate and label two variables
drawnorm sat motivation, ///
n(250) ///
means(.007, -.053) ///
sds(1.01, 1.02) ///
corr(1, .035, 1) cstorage(lower) ///
clear seed(1)
label var sat "SAT"
label var motivation "Motivation"
// Only the 15 per cent at the top are admitted
gen admission_sc = sat + motivation
_pctile admission_sc, percentiles(85)
gen admission = (admission_sc > r(r1))
label var admission "Admission status"
label define admission 1 "Admitted applicant" ///
0 "Rejected applicant"
label val admission admission
drop admission_sc
// Plot as in Morgan and Winship, p. 67:
twoway (scatter motivation sat if admission == 1) ///
(scatter motivation sat if admission == 0) ///
, ///
legend(label(1 "Admitted applicants") ///
label(2 "Rejected applicants") ///
pos(5) ring(0)) ///
ylabel(none) xlabel(none) ///
name(collider1, replace)
// Enhanced plot with fitted lines and correlations
quietly cor motivation sat
local r_overall = round(r(rho), .01)
quietly cor motivation sat if admission == 1
local r_admitted = round(r(rho), .01)
quietly cor motivation sat if admission == 0
local r_rejected = round(r(rho), .01)
twoway (scatter motivation sat if admission == 1) ///
(scatter motivation sat if admission == 0) ///
(lfit motivation sat) ///
(lfit motivation sat if admission == 1) ///
(lfit motivation sat if admission == 0) ///
, ///
legend(label(1 "Admitted applicants") ///
label(2 "Rejected applicants") ///
label(3 "Overall fit, {it:r} = `r_overall'") ///
label(4 "Fit for admitted, {it:r} = `r_admitted'") ///
label(5 "Fit for rejected, {it:r} = `r_rejected'") ///
pos(5) ring(0)) ///
ylabel(none) xlabel(none) ytitle("Motivation")
What the Figures show is that the very small correlation between motivation and SAT score for the overall group turns out to be much larger when conditioning for admission status.
This also shows in an OLS regression:
regress motivation sat, beta
estimates store m1
regress motivation sat admission, beta
estimates store m2
estimates table m1 m2, b(%7.2f) se(%7.2f) stats(N) label
----------------------------------------------
Variable | m1 m2
-------------------------+--------------------
SAT | 0.10 -0.21
| 0.06 0.06
Admission status | 1.60
| 0.17
Constant | -0.12 -0.36
| 0.06 0.06
-------------------------+--------------------
N | 250 250
----------------------------------------------
legend: b/se
Cole et al. (2010) present additional illustrations for this problem.References
Cole, Stephen R., Robert W. Platt, Enrique F. Schisterman, Haitao Chu, Daniel Westreich, David Richardson, and Charles Poole. 2010. "Illustrating Bias Due to Conditioning on a Collider." International Journal of Epidemiology 39(2):417-420. doi: 10.1093/ije/dyp334
Morgan, Stephen L., and Christopher Winship. 2007. Counterfactuals and Causal Inference. Methods and Principles for Social Research. Cambridge University Press.
Pearl, Judea. 2009. Causality. Models, Reasoning, and Inference, 2nd ed. Cambridge University Press.
Labels:
_pctile,
drawnorm,
estimates table,
Legends,
local,
round(),
Simulation,
Textbooks,
twoway lfit,
twoway scatter
Oct 5, 2012
Random graphs (4): Plotting intercept and slope variation
use V3 V51 ISCO88 using "C:\Users\User\work\data sets\issp - work orientations\work orientations iii (2005)\ZA4350_F1.dta", clear // Generate variables gen jobsatisf = 7 - V51 iskoisei isei, isko(ISCO88) drop V51 ISCO88 preserve // Calculating intercept and slope variation using OLS statsby inter = _b[_cons] /// slope = _b[isei] /// , by(V3) /// saving(ols, replace): /// regress jobsatisf isei merge m:1 V3 using ols drop _merge // Visualizing intercept and slope variation gen yhat_ols = inter + slope*isei separate jobsatisf, by(V3) separate yhat_ols, by(V3) twoway (line yhat_ols1-yhat_ols43 isei, sort(V3 isei)) /// (lfit jobsatisf isei, clwidth(vvthick) clcolor(black)) /// , legend(off) ytitle("Job satisfaction") xtitle("ISEI") /// xlabel(16 25 50 75 90) ylabel(,format(%6.1f)) /// caption("{it:Source:} ISSP 2005 (Work Orientations III), own calculations", span) /// name(one, replace) restore // Calculating intercept and slope variation using a random effects model center isei, inplace mixed jobsatisf isei || V3: isei, var predict u1 u0, reffects // Visualizing intercept and slope variation gen predRandomSlope = (_b[_cons] + u0) + ((_b[isei] + u1) * isei) twoway (line predRandomSlope isei, connect(ascending) sort(V3 isei)), /// ytitle("Job satisfaction") /// xtitle("ISEI (centered)") /// xlabel(-25 0 25 50) /// ylabel(,format(%6.1f)) /// caption("{it:Source:} ISSP 2005 (Work Orientations III), own calculations", span) /// name(two, replace)
Subscribe to:
Posts (Atom)









