Effective graphs for data display

Introduction
Figure 1
Figure 2
Figure 3
Figure 4
Figure 5
Figure 6
Figure 7
Figure 8
Figure 9
Figure 10
Figure 11
Figure 12

Introduction

This page contains Stata code and output, equivalent to the examples in the paper. In all cases, the data are random samples from the freely-available NHANES dataset, and the variables used are;

BPXSAR: systolic blood pressure (mmHg)
BPXDAR: diastolic blood pressure (mmHg)
BPXDI1, BPXDI2: two diastolic blood pressure readings
race_ethc: race/ethnicity, coded as Hispanic, White non-Hispanic, Black non-Hispanic and Other
DR1TFOLA: folate intake (\(\mu\)g/day)
RIAGENDR: sex, coded as Male/Female
BMXBMI: body mass index (kg/m\(^2\))
RIDAGEY: age (years)

To run the code, the following packages must first be installed;

vioplot: ssc install vioplot
stripplot: ssc install stripplot
spineplot: ssc install spineplot

Following this, cut-and-paste the code for each figure into an R session.

Figure 1

insheet using http://faculty.washington.edu/kenrice/heartgraphs/nhanessmall.csv, clear
label variable bpxsar "Systolic Blood Pressure (mmHg)"
stripplot bpxsar

Figure 1a

stripplot bpxsar, stack center

Figure 1b

Figure 2

insheet using http://faculty.washington.edu/kenrice/heartgraphs/nhanesmedium.csv, clear
label variable bpxsar "Systolic Blood Pressure (mmHg)"
stripplot bpxsar

Figure 2a

stripplot bpxsar, stack center

Figure 2b

stripplot bpxsar, width(1) stack center

Figure 2c

Figure 3

insheet using http://faculty.washington.edu/kenrice/heartgraphs/nhaneslarge.csv, clear
label variable bpxsar "Systolic Blood Pressure (mmHg)"
vioplot bpxsar, horizontal

Figure 3a

histogram bpxsar

Figure 3b

Figure 4

insheet using http://faculty.washington.edu/kenrice/heartgraphs/nhaneslarge.csv, clear
generate normotens = bpxsar <= 140
label define normotensl 1 "<= 140 mmHg" 0 "> 140 mmHg" 
label values normotens normotensl
generate sample = 1
graph hbar (percent) sample, over(normotens) asyvars stack

Figure 4a

insheet using http://faculty.washington.edu/kenrice/heartgraphs/nhaneslarge.csv, clear
generate normotens = bpxsar <= 140
collapse (mean) est=normotens (sebinomial) stderr=normotens
generate loCI = est - 1.96*stderr
generate hiCI = est + 1.96*stderr
generate x = 1
label variable x "Estimate"
label variable loCI "95% Conf Int"
label variable hiCI "95% Conf Int"
twoway (scatter x est) (rspike loCI hiCI x, horizontal), ///
xscale( r(0 1) ) xtick( 0(0.2)1 ) xlabel( 0(0.2)1) yscale(off)

Figure 4b

insheet using http://faculty.washington.edu/kenrice/heartgraphs/nhaneslarge.csv, clear
generate sample = 1
graph hbar (percent) sample, over(race_ethc) asyvars stack

Figure 4c

insheet using http://faculty.washington.edu/kenrice/heartgraphs/nhaneslarge.csv, clear
generate x=1
collapse (count) x, by(race_ethc)
generate est = x/1000
generate stderr = sqrt(est*(1-est)/1000)
generate loCI = est - 1.96*stderr
generate hiCI = est + 1.96*stderr
generate study = _n
label variable study "Estimate"
label variable loCI "95% Conf Int"
label variable hiCI "95% Conf Int"
twoway (scatter study est) (rspike loCI hiCI study, horizontal), ///
xscale( r(0 1) ) xtick( 0(0.2)1 ) xlabel( 0(0.2)1) /// 
ylabel( 1 "Black non-Hispanic" 2 "Hispanic" 3 "Other" 4 "White non-Hispanic", angle(0)) ///
ytitle("")

Figure 4d

Figure 5

insheet using http://faculty.washington.edu/kenrice/heartgraphs/nhaneslarge.csv, clear
generate normotens = bpxsar <= 140
label define normotensl 1 "<= 140 mmHg" 0 "> 140 mmHg" 
label values normotens normotensl
generate x = 1
graph pie x, over(normotens)

Figure 5a

insheet using http://faculty.washington.edu/kenrice/heartgraphs/nhaneslarge.csv, clear
generate x = 1
graph pie x, over(race_ethc)

Figure 5b

Figure 6

insheet using http://faculty.washington.edu/kenrice/heartgraphs/nhanessmall.csv, clear
label variable dr1tfola "Folate intake (ug/day)"
label define sexl 1 "Male" 2 "Female" 
label values riagendr sexl
stripplot dr1tfola, vertical over(riagendr)

Figure 6a

insheet using http://faculty.washington.edu/kenrice/heartgraphs/nhanesmedium.csv, clear
label variable dr1tfola "Folate intake (ug/day)"
label define sexl 1 "Male" 2 "Female" 
label values riagendr sexl
stripplot dr1tfola, width(20) vertical stack center over(riagendr)

Figure 6b

insheet using http://faculty.washington.edu/kenrice/heartgraphs/nhaneslarge.csv, clear
label variable dr1tfola "Folate intake (ug/day)"
label define sexl 1 "Male" 2 "Female" 
label values riagendr sexl
vioplot dr1tfola, vertical over(riagendr) ytitle("Folate intake (ug/day)")

Figure 6c

Figure 7

insheet using http://faculty.washington.edu/kenrice/heartgraphs/nhanesmedium.csv, clear
label variable bpxsar "Systolic Blood Pressure (mmHg)"
label variable ridageyr "Age (years)"
twoway (scatter bpxsar ridageyr, msymbol(Oh)) (lfitci bpxsar ridageyr), legend(off) ytitle("Systolic Blood Pressure (mmHg)")

Figure 7a

insheet using http://faculty.washington.edu/kenrice/heartgraphs/nhaneslarge.csv, clear
label variable bpxsar "Systolic Blood Pressure (mmHg)"
label variable ridageyr "Age (years)"
mkspline spage = ridageyr, cubic
quietly regress bpxsar spage*
predict fit, xb
predict stderr, stdp
generate lb = fit - 1.96*stderr
generate ub = fit + 1.96*stderr
sort ridageyr 
twoway scatter bpxsar ridageyr || rarea lb ub ridageyr || line fit ridageyr, legend(off) ytitle("Systolic Blood Pressure (mmHg)")

Figure 8

insheet using http://faculty.washington.edu/kenrice/heartgraphs/nhaneslarge.csv, clear
egen sbpcut = cut(bpxsar), at(0 110 140 250) label
encode race_ethc, generate(race_ethc2)
label variable sbpcut "Category of Sys Blood Pressure"
label variable race_ethc2 "Race/Ethnicity"
spineplot sbpcut race_ethc2, xlabel(, angle(v) axis(2))

Figure 8a

insheet using http://faculty.washington.edu/kenrice/heartgraphs/nhaneslarge.csv, clear
generate x=1
egen sbpcut = cut(bpxsar), at(0 110 140 250) label
collapse (count) x, by(race_ethc sbpcut)
egen prop = pc(x), by(race_ethc)
replace prop = prop/100
encode race_ethc, generate(race_ethc2)
label variable prop "Proportion"
label variable race_ethc2 "Race/Ethnicity"
stripplot prop, over(race_ethc2) separate(sbpcut) vertical xlabel(, angle(v)) legend( title("Category of Sys Blood Pressure"))

Figure 8b

Figure 9

insheet using http://faculty.washington.edu/kenrice/heartgraphs/nhaneslarge.csv, clear
egen dbpcut1 = cut(bpxdi1), at(0 60 70 80 120) label
egen dbpcut2 = cut(bpxdi2), at(0 60 70 80 120) label
generate x=1
collapse (count) x, by(dbpcut1 dbpcut2) 
twoway (scatter dbpcut1 dbpcut2 [weight=x], msize(1.5) msymbol(S)) ///
(function y=x, range(-1 4) ), ///
aspectratio(1) ///
xlabel(0 "0-" 1 "60-" 2 "70-" 3 "80-") ///
ylabel(0 "0-" 1 "60-" 2 "70-" 3 "80-") /// 
xscale( range(-1 4)) ///
yscale( range(-1 4)) ///
xtitle("Category of Dias BP, Measure #2") ///
ytitle("Category of Dias BP, Measure #1") ///
legend(off)

Figure 10

insheet using http://faculty.washington.edu/kenrice/heartgraphs/nhanesmedium.csv, clear
generate htn = bpxdar > 70
label variable ridageyr "Age (years)"
label variable htn "Diastolic BP (mmHg)/Proportion"
label define htnl 1 "> 70 mmHg" 0 "<= 70 mmHg" 
label values htn htnl
logistic htn ridageyr
predict fitted, xb
predict fitted_se, stdp
generate est = exp(fitted)/(1+exp(fitted))
generate hici = exp(fitted+fitted_se*1.96)/(1+exp(fitted+fitted_se*1.96))
generate loci = exp(fitted-fitted_se*1.96)/(1+exp(fitted-fitted_se*1.96))
sort ridageyr
stripplot ridageyr, stack center over(htn) 
graph addplot rarea hici loci ridageyr
graph addplot line est ridageyr, yscale( range(0 1)) ylabel(0 "<= 70mmHg 0" 0.5 1 ">70mmHg 1")

Figure 10a

insheet using http://faculty.washington.edu/kenrice/heartgraphs/nhanesmedium.csv, clear
generate htn = bpxdar > 70
label variable ridageyr "Age (years)"
label variable htn "Dichotomized Diastolic BP (mmHg), Proportion"
label define htnl 1 "> 70 mmHg" 0 "<= 70 mmHg" 
label values htn htnl
mkspline spage = ridageyr, cubic
quietly logistic htn spage*
predict fitted, xb
predict fitted_se, stdp
generate est = exp(fitted)/(1+exp(fitted))
generate hici = exp(fitted+fitted_se*1.96)/(1+exp(fitted+fitted_se*1.96))
generate loci = exp(fitted-fitted_se*1.96)/(1+exp(fitted-fitted_se*1.96))
sort ridageyr
stripplot ridageyr, stack center over(htn) 
graph addplot rarea hici loci ridageyr
graph addplot line est ridageyr, yscale( range(0 1)) ylabel(0 "<= 70mmHg 0" 0.5 1 ">70mmHg 1")

Figure 10b

Figure 11

insheet using http://faculty.washington.edu/kenrice/heartgraphs/nhaneslarge.csv, clear
egen dbpcut = cut(bpxdar), at(0 60 70 80 120) label
forvalues i = 0(1)3{
    generate level`i' = dbpcut==`i'
    quietly logit level`i' ridageyr
    predict fitxb`i', xb
    generate fitp`i' = exp(fitxb`i')/(1+exp(fitxb`i'))
    local temp : label (dbpcut) `i'
    label variable fitp`i' "`temp'"
    }
    
sort ridageyr
label variable ridageyr "Age (years)"
line fitp* ridageyr, ytitle("Proportion, based on log'c regression")

Figure 11a

insheet using http://faculty.washington.edu/kenrice/heartgraphs/nhaneslarge.csv, clear
egen dbpcut = cut(bpxdar), at(0 60 70 80 120) label
forvalues i = 0(1)2{
    generate level`i' = dbpcut>`i'
    quietly logit level`i' ridageyr
    predict fitxb`i', xb
    generate fitp`i' = exp(fitxb`i')/(1+exp(fitxb`i'))
    }
    
label variable ridageyr "Age (years)"
label variable fitp0 "> 60 mmHg"
label variable fitp1 "> 70 mmHg"
label variable fitp2 "> 80 mmHg"
sort ridageyr
line fitp* ridageyr, ytitle("Proportion, based on log'c regression")

Figure 11b

Figure 12

insheet using http://faculty.washington.edu/kenrice/heartgraphs/nhaneslarge.csv, clear
egen agecut = cut(ridageyr), at(0 30 55 100) label
label variable bmxbmi "BMI (kg/m^2)"

twoway (scatter bpxsar bmxbmi, msymbol(o) mcolor(gs14) subtitle("Age Category:", prefix)) ///
(lowess bpxsar bmxbmi if riagendr==1, lwidth(thick)) ///
(lowess bpxsar bmxbmi if riagendr==2, lwidth(thick) lpattern(dash)) ///
, by(agecut, cols(3) compact note("")) ///
ytitle("Systolic Blood Pressure (mmHg)") legend(order(2 "Male" 3 "Female"))