clear // Install rcall command *github install haghish/rcall, stable // Clear memory of R session rcall clear //#install.packages("essurvey"); rcall: library(essurvey) // Download data rcall: download_rounds(c(9), /// ess_email = "REGISTERED E-MAIL ADDRESS HERE", /// output_dir = "data", /// format = 'stata') // Open data set local vars cntry gndr polintr dweight use `vars' using "data/ESS9/ESS9e01_2.dta", clear // Variables generate interest = (inlist(polintr, 1, 2)) replace interest = . if missing(polintr) generate female = (gndr == 2) // Calculate gender gap per country tempname foo postfile `foo' str50 cntry men women gap gap_se using deleteme.dta, replace levelsof cntry, local(cntry) foreach i of local cntry { qui regress interest i.female [pw = dweight] if cntry == `"`i'"' qui margins i.female, post local men = _b[0.female] * 100 local women = _b[1.female] * 100 qui regress interest i.female [pw = dweight] if cntry == `"`i'"' qui margins r.female, post local gap = _b[r1vs0.female] * -100 local gap_se = _se[r1vs0.female] * -100 post `foo' ("`i'") (`men') (`women') (`gap') (`gap_se') } postclose `foo' // Open data set with gender gap per country use deleteme.dta, clear kountry cntry, from(iso2c) egen country = rank(gap), unique labmask country, value(NAMES_STD) generate ub = gap + 1.96 * gap_se generate lb = gap - 1.96 * gap_se generate y1 = -women generate y2 = 0 generate y3 = -women + men generate x1 = country - 0.2 generate x2 = country + 0.2 mylabels 0(20)60, myscale(-@) local(show) local spaces = 30 * " " twoway (scatter y1 country, msymbol(p) yaxis(1 2)) /// (rbar y1 y2 x1, blcolor(gs10%40) bfcolor(gs10%40) barwidth(0.4)) /// (rbar y1 y3 x2, blcolor(gs5%40) bfcolor(gs5%40) barwidth(0.4)) /// (rbar gap y2 x2, blcolor(gs3%40) bfcolor(gs3%40) barwidth(0.4)) /// (rspike ub lb x2) /// , ylabel(`show', grid ang(h)) /// ylabel(0(10)40, grid ang(h) axis(2)) /// ytitle("% interested in politics") /// ytitle("`spaces' Gender gap", axis(2)) /// xlabel(1/19, noticks ang(45) valuelabel labsize(small)) /// legend(order(2 "Men" 3 "Women" 4 "Gender gap" 5 "Gender gap 95% CI") /// col(4) pos(12) ring(0)) name(gendergap, replace)
Showing posts with label Confidence intervals. Show all posts
Showing posts with label Confidence intervals. Show all posts
May 19, 2020
Random graphs (146): Bar graph
Sep 29, 2015
Scaling predictors for models with interactions
The scaling of predictors is an important issue for the interpretation of models with interaction terms. The Table below presents a simple OLS regression model of happiness regressed on age, sex, and their interaction.
When none of the predictors are standardized and 0 is not a meaningful value on a predictor, interpretation of the regression coefficients can become awkward. For instance, in the model with none of the predictors standardized, the coefficient for female – .348 – states that women score .348 points higher in terms of happiness than men among those who are 0 years old. The coefficient for age – -.004 – indicates that with every additional year, men's happiness decreases by .004 points. The coefficient for the interaction term – -.008 – denotes that women's happiness decrease more than men's with each additional year of age. The intercept is the average happiness of men who are 0 years old.
A first step towards giving the coefficients from an interaction model an easier interpretation is to standardize continuous variables. For the second model shown in the Table, age was standardized by subtracting its mean and dividing it by its standard deviation. Thus, for the standardized variable, the value of 0 refers to the average age, and a change of 1 refers to a change of one standard deviation. The age coefficient (-.07) now indicates that for every additional standard deviation of age, men's happiness decreases by .07. For men and women of average age, the happiness differential is .02. Women's happiness decreases by .15 per standard deviation of age stronger than men's. The intercept shows the average happiness of average-age men.
While it is a rather controversial technique, it is in principle also possible to standardize the dummy variable indicating the difference between men and women. Then, the intercept is the overall mean of happiness. However, an increase of 1 in the gender coefficient becomes difficult to interpret, as there is no such thing as a one-standard deviation increase in gender. When the dummy variable is approximately evenly distributed, then the coefficient reflects half of the difference in the outcome between the two groups.
Gelman and Hill (2007, p. 56) suggest to standardize variables by two rather than one standard deviation to keep binary and continuous variables on a roughly common scale. An increase of one in age then refers to the difference between one standard deviation below the mean and one standard deviation above the mean; an increase of one in age reflects roughly the difference between men and women.
As the plot above shows, none of these different scaling decisions affect the form of the interaction.
When none of the predictors are standardized and 0 is not a meaningful value on a predictor, interpretation of the regression coefficients can become awkward. For instance, in the model with none of the predictors standardized, the coefficient for female – .348 – states that women score .348 points higher in terms of happiness than men among those who are 0 years old. The coefficient for age – -.004 – indicates that with every additional year, men's happiness decreases by .004 points. The coefficient for the interaction term – -.008 – denotes that women's happiness decrease more than men's with each additional year of age. The intercept is the average happiness of men who are 0 years old.
A first step towards giving the coefficients from an interaction model an easier interpretation is to standardize continuous variables. For the second model shown in the Table, age was standardized by subtracting its mean and dividing it by its standard deviation. Thus, for the standardized variable, the value of 0 refers to the average age, and a change of 1 refers to a change of one standard deviation. The age coefficient (-.07) now indicates that for every additional standard deviation of age, men's happiness decreases by .07. For men and women of average age, the happiness differential is .02. Women's happiness decreases by .15 per standard deviation of age stronger than men's. The intercept shows the average happiness of average-age men.
While it is a rather controversial technique, it is in principle also possible to standardize the dummy variable indicating the difference between men and women. Then, the intercept is the overall mean of happiness. However, an increase of 1 in the gender coefficient becomes difficult to interpret, as there is no such thing as a one-standard deviation increase in gender. When the dummy variable is approximately evenly distributed, then the coefficient reflects half of the difference in the outcome between the two groups.
Gelman and Hill (2007, p. 56) suggest to standardize variables by two rather than one standard deviation to keep binary and continuous variables on a roughly common scale. An increase of one in age then refers to the difference between one standard deviation below the mean and one standard deviation above the mean; an increase of one in age reflects roughly the difference between men and women.
Reference
Gelman, Andrew, and Jennifer Hill. 2007. Data Analysis Using Regression and Multilevel/Hierarchical Models. Cambridge University Press.use agea gndr happy using ESS1e06.3_F1.dta, clear
// Set up variables
rename agea age
recode gndr (1 = 0) (2 = 1) (.a = .), gen(female)
label var female "Female sex (ref. male)"
label var age "Age"
label var happy "Happiness"
// Listwise deletion
drop if missing(female, age, happy)
// Standardize variables
qui sum female
gen female_1z = (female - r(mean)) / r(sd)
gen female_2z = (female - r(mean)) / (2*r(sd))
qui sum age
gen age_1z = (age - r(mean)) / r(sd)
gen age_2z = (age - r(mean)) / (2*r(sd))
// Fit models
reg happy c.age##female
estimates store unstandardized
reg happy c.age_1z##female
estimates store cont_standardized
reg happy c.age_1z##c.female_1z
estimates store standardized_1z
reg happy c.age_2z##c.female_2z
estimates store standardized_2z
// Create table
esttab unstandardized cont_standardized standardized_1z standardized_2z, ///
rename(age_1z age ///
age_2z age ///
female_1z 1.female ///
female_2z 1.female ///
c.age_1z#c.female_1z 1.female#c.age ///
c.age_2z#c.female_2z 1.female#c.age ///
1.female#c.age_1z 1.female#c.age) ///
drop(0b.female 0b.female#co.age 0b.female#co.age_1z) ///
coeflabel(age "Age" ///
1.female "Female (ref. male)" ///
1.female#c.age "Age X Female" ///
_cons "Intercept") ///
mtitles("Unstandardized" "Age standardized (1 SD)" ///
"Standardized (1 SD)" "Standardized (2 SD)") ///
nonumbers varwidth(18) modelwidth(24) ///
se r2 obslast
// Plot models
estimates restore unstandardized
qui sum age
local age_minsd = r(mean) - r(sd)
local age_mean = r(mean)
local age_plusd = r(mean) + r(sd)
qui margins, at(c.age = (`age_minsd' `age_mean' `age_plusd') ///
female = (0 1)) vsquish
qui marginsplot, recastci(rarea) ciopts(color(gs12)) ///
xlabel(`age_minsd' "-1 SD" `age_mean' "Average age" `age_plusd' "+1 SD") ///
title("Unstandardized variables", size(*.9)) ///
ytitle("Predicted happiness") ///
xtitle("") ///
plotopts(msymbol(none)) /// // Turn off markers
plot1opts(lpattern(dash)) /// // Define line types here
plot2opts(lpattern(solid)) ///
legend(subtitle(Gender, size(*.8)) ///
order(3 "Males" 4 "Females") ring(0) pos(7) ///
size(*.8) symysize(*.7) symxsize(*.7)) ///
name(unstandardized, replace)
estimates restore cont_standardized
qui margins, at(c.age_1z = (-1 0 1) ///
female = (0 1)) vsquish
qui marginsplot, recastci(rarea) ciopts(color(gs12)) ///
xlabel(-1 "-1 SD" 0 "Average age" 1 "+1 SD") ///
title("Age standardized at 1 SD", size(*.9)) ///
ytitle("Predicted happiness") ///
xtitle("") ///
plotopts(msymbol(none)) /// // Turn off markers
plot1opts(lpattern(dash)) /// // Define line types here
plot2opts(lpattern(solid)) ///
legend(subtitle(Gender, size(*.8)) ///
order(3 "Males" 4 "Females") ring(0) pos(7) ///
size(*.8) symysize(*.7) symxsize(*.7)) ///
name(cont_standardized, replace)
estimates restore standardized_1z
qui sum female
local fem_hi = (1 - r(mean)) / r(sd)
local fem_lo = (0 - r(mean)) / r(sd)
qui margins, at(c.age_1z = (-1 0 1) ///
c.female_1z = (`fem_lo' `fem_hi')) vsquish
qui marginsplot, recastci(rarea) ciopts(color(gs12)) ///
xlabel(-1 "-1 SD" 0 "Average age" 1 "+1 SD") ///
title("All variables standardized at 1 SD", size(*.9)) ///
ytitle("Predicted happiness") ///
xtitle("") ///
plotopts(msymbol(none)) /// // Turn off markers
plot1opts(lpattern(dash)) /// // Define line types here
plot2opts(lpattern(solid)) ///
legend(subtitle(Gender, size(*.8)) ///
order(3 "Males" 4 "Females") ring(0) pos(7) ///
size(*.8) symysize(*.7) symxsize(*.7)) ///
name(standardized_1z, replace)
estimates restore standardized_2z
qui sum female
local fem_hi = (1 - r(mean)) / (2*r(sd))
local fem_lo = (0 - r(mean)) / (2*r(sd))
qui margins, at(c.age_2z = (-.5 0 .5) ///
c.female_2z = (`fem_lo' `fem_hi')) vsquish
qui marginsplot, recastci(rarea) ciopts(color(gs12)) ///
xlabel(-.5 "-1 SD" 0 "Average age" .5 "+1 SD") ///
title("All variables standardized at 2 SD", size(*.9)) ///
ytitle("Predicted happiness") ///
xtitle("") ///
plotopts(msymbol(none)) /// // Turn off markers
plot1opts(lpattern(dash)) /// // Define line types here
plot2opts(lpattern(solid)) ///
legend(subtitle(Gender, size(*.8)) ///
order(3 "Males" 4 "Females") ring(0) pos(7) ///
size(*.8) symysize(*.7) symxsize(*.7)) ///
name(standardized_2z, replace)
graph combine unstandardized cont_standardized standardized_1z standardized_2z, ///
row(2) col(2) ysize(7) name(combined, replace)
May 24, 2015
Random graphs (48): Interaction plot with overlaid data points
// Simulate data
clear
set seed 1
set obs 500
generate e = 0 + (500 - 0) * runiform() // To generate random variates over the
generate x1 = 0 + (800 - 0) * runiform() // interval [a,b), a+(b-a)*runiform()
generate x2 = round(0 + (1 - 0) * runiform()) // Binary variable
generate y = (x1 + x2 + (x1*x2) + e) / 1000
// Calculations for Aiken & West (1991) style plot
regress y c.x1##x2
qui sum x1
local x1_minsd = r(mean) - (2*r(sd))
local x1_mean = r(mean)
local x1_plusd = r(mean) + (2*r(sd))
// Calculate predicted values for plot
margins, at(c.x1 = (`x1_minsd' `x1_mean' `x1_plusd') ///
c.x2 = (0 1)) vsquish
// Plot
qui marginsplot, recastci(rarea) ciopts(color(gs10)) ///
title("Interaction plot with overlaid data points") ///
ytitle("y") ///
ylabel(, format(%6.1f)) ///
xtitle("") ///
plotopts(msymbol(none)) /// // Turn off markers
plot1opts(lpattern(longdash)) /// // Define line types here
plot2opts(lpattern(solid)) ///
addplot(scatter y x1 ///
, symbol(o) ///
xlabel(`x1_minsd' "-2 SD" /// // The addplot seems to override
`x1_mean' "Average x1" /// // the regular axis label
`x1_plusd' "+2 SD") /// // command
legend(subtitle(x2) /// // And the legend command
order(3 "x2 = 0" 4 "x2 = 1" 2 "95 % CI"))) ///
legend(pos(5) ring(0)) /// // But not completely
name(plot, replace)
Apr 20, 2015
Random graphs (45): Re-creating -graph dot- using -twoway dot-
This plot recreates this one, but adds confidence intervals.
preserve
// Calculate country-specific means to be plotted
tempname foo
tempname idealage
postfile `foo' quintile incocomp incocomp_lb incocomp_ub using `foo2', replace
levelsof(quintile), local(quintile)
foreach x of local quintile {
qui reg incocomp if quintile == `x', cluster(cntry)
local incocomp = _b[_cons]
local incocomplb = _b[_cons] - (1.96 * _se[_cons])
local incocompub = _b[_cons] + (1.96 * _se[_cons])
post `foo' (`x') (`incocomp') (`incocomplb') (`incocompub')
}
postclose `foo'
// Plot country-specific means
use `foo2', clear
twoway (dot incocomp quintile) /// -twoway scatter- doesn't have the lines
(rcap incocomp_lb incocomp_ub quintile) ///
, yscale(range(1.5 3.0)) ///
ylabel(1.5 (.25) 3.0, format(%6.2f)) ///
legend(off)
title("Income quintiles")
xlabel(1 `""1" "(poorest)""' /// Double-compound quotes allow for the line break
2 3 4 5 `""5" "(richest)"') ///
fxsize(40) /// Forces x-size to be 40 percent of its original size
xtitle("") ///
xscale(range(0.2 5.8)) /// Making this a bit bigger makes room for the labels
name(byquintile_hor, replace)
restore
// Calculate country-specific means to be plotted 2
preserve
tempname foo
tempname idealage
postfile `foo' str2 cntry incocomp incocomp_lb incocomp_ub using `foo2', replace
levelsof(cntry), local(cntry)
foreach x of local cntry {
qui reg incocomp if cntry == "`x'"
local incocomp = _b[_cons]
local incocomplb = _b[_cons] - (1.96 * _se[_cons])
local incocompub = _b[_cons] + (1.96 * _se[_cons])
post `foo' ("`x'") (`incocomp') (`incocomplb') (`incocompub')
}
postclose `foo'
use `foo2', clear
// Create sorted-by-size variable
egen order_ = rank(-incocomp), unique
labmask order_, value(cntry)
// Plot country-specific means 2
twoway (dot incocomp order_) ///
(rcap incocomp_lb incocomp_ub order_) ///
, yscale(off) /// Not necessary, as we're using the one from the other graph
xlabel(1/23, val alt) ///
title("Countries") ///
xtitle("") ///
legend(off) ///
/* graphregion(margin(b=16)) */ /// Could help as -ycommon- isn't accounting for the different label sizes
name(bycntry_hor, replace)
// Combine plots
graph combine byquintile_hor bycntry_hor, ///
l1title("Mean income comparison orientation") ///
ycommon row(1) ///
imargin(zero) // Sets margin between both plots to 0
restore
Nov 30, 2014
Random graphs (37): Proportions with confidence intervals
use acountry asex aage ayear a602 a612 aweight using GGS_Wave1_V.4.2.dta, clear // Create outcome variable recode a612 (1 2 1601 = 1 "(Likely) infertile") /// (3 4 1501 1602 = 0 "(Likely) fertile") /// (2001 . .a .b .c = .x "Don't know"), /// gen(infertility) replace infertility = 0 if a602 == 1 // Also fertile when currently pregnant // Create country identifier gen cntry = "" replace cntry = "BG" if acountry == 11 replace cntry = "RU" if acountry == 12 replace cntry = "GE" if acountry == 13 replace cntry = "DE" if acountry == 14 replace cntry = "FR" if acountry == 15 replace cntry = "HU" if acountry == 16 replace cntry = "IT" if acountry == 17 replace cntry = "NL" if acountry == 18 replace cntry = "RO" if acountry == 19 replace cntry = "NO" if acountry == 20 replace cntry = "AT" if acountry == 21 replace cntry = "EE" if acountry == 22 replace cntry = "BE" if acountry == 23 replace cntry = "AU" if acountry == 24 replace cntry = "LT" if acountry == 25 replace cntry = "PL" if acountry == 26 replace cntry = "CZ" if acountry == 28 // Recode gender variable gen sex = "" replace sex = "Men" if asex == 1 replace sex = "Women" if asex == 2 // Drop unnecessary countries drop if inlist(cntry, "IT", "AU", "NL") // Prevalence plot by country preserve // Drop unnecesary age groups keep if aage <= 39 keep if aage >= 34 svyset [pweight = aweight] // Declare weight // Create temporary files and postfile tempname foo tempname infert postfile `foo' str2 cntry str5 sex perc_infert perc_infert_ll perc_infert_ul using `infert' , replace levelsof cntry, local(levels1) levelsof sex, local(sex1) // Calculate proportions by country and sex foreach sex of local sex1 { foreach country of local levels1 { // Calculate proportions capture proportion infertility if cntry == "`country'" & sex == "`sex'" *matrix list r(table) matrix coefs = r(table) local perc = coefs[1,2] * 100 local perc_ll = coefs[5,2] * 100 local perc_ul = coefs[6,2] * 100 di "`country'" _skip(2) `perc_ll' _skip(2) `perc' _skip(2) `perc_ul' if "`perc'" != "" { // Make sure that the loop doesn't break if empty post `foo' ("`country'") ("`sex'") (`perc') (`perc_ll') (`perc_ul') } } } postclose `foo' // Use posted data set use `infert', clear // Sort countries by size for women egen order_ = rank(-perc_infert) if sex == "Women", unique labmask order_, value(cntry) // Plot for women twoway (dot perc_infert order_ if sex == "Women", horizontal) /// (rcap perc_infert_ll perc_infert_ul order_ if sex == "Women", horizontal) /// , ylabel(1/14, val) legend(off) /// ytitle("") xtitle("% of {bf:women} aged 35{c 0150}39 reporting (suspected) infertility") /// /*caption("{it:Note:} Error bars denote 95 % confidence intervals", size(small) span)*/ /// name(infert3539women, replace) // Drop women's order drop order_ // Sort countries by size for men egen order_ = rank(-perc_infert) if sex == "Men", unique labmask order_, value(cntry) // Plot for men twoway (dot perc_infert order_ if sex == "Men", horizontal) /// (rcap perc_infert_ll perc_infert_ul order_ if sex == "Men", horizontal) /// , ylabel(1/13, val) legend(off) /// ytitle("") xtitle("% of {bf:men} aged 35{c 0150}39 reporting (suspected) infertility") /// /*caption("{it:Note:} Error bars denote 95 % confidence intervals", size(small) span)*/ /// name(infert3539men, replace) // Combine graphs graph combine infert3539women infert3539men /// , ysize(8) col(1) xcommon /// caption("{it:Source:} Generations and Gender Survey Wave 1 (2005{c 0150}2011}" /// "{it:Note:} Error bars denote 95 % confidence intervals", size(vsmall) span) restore// Plot infertility by age // Recode age into groups recode aage (20/24 = 1 "20–24") /// (25/29 = 2 "25–29") /// (30/34 = 3 "30–34") /// (35/39 = 4 "35–39") /// (40/45 = 5 "40–45") /// (1/19 = .a "Below 20 y") /// (46/85 = .b "46+ y") /// , gen(age) // Save valuelabels in temporary do-file tempfile valuelabels label save age using `valuelabels' svyset [pweight = aweight] // Declare weight // Create temporary files and postfile tempname foo tempname infertage postfile `foo' str2 cntry str5 sex age perc_infert perc_infert_ll perc_infert_ul using `infertage' , replace levelsof cntry, local(levels1) levelsof sex, local(sex1) // Calculate proportion by age group, country, and sex foreach sex of local sex1 { foreach country of local levels1 { forvalues x = 1/5 { // Calculate proportions capture proportion infertility if age == `x' & cntry == "`country'" & sex == "`sex'" *matrix list r(table) matrix coefs = r(table) local perc = coefs[1,2] * 100 local perc_ll = coefs[5,2] * 100 local perc_ul = coefs[6,2] * 100 di "`country'" _skip(2) `x' _skip(2) `perc_ll' _skip(2) `perc' _skip(2) `perc_ul' if "`perc'" != "" { // Make sure that the loop doesn't break if post `foo' ("`country'") ("`sex'") (`x') (`perc') (`perc_ll') (`perc_ul') } } } } postclose `foo' // Use posted values use `infertage', clear drop if perc_infert == . // Drop empty rows // Define and set saved value labels again do `valuelabels' label val age age // Plot for women by country and age twoway (rarea perc_infert_ll perc_infert_ul age if sex == "Women", fcolor(gs15) cmissing(no)) /// (line perc_infert age if sex == "Women", cmissing(no)) /// , by(cntry, cols(3) legend(off) note("") /// caption("{it:Source:} Generations and Gender Survey Wave 1 (2005{c 0150}2011)" /// "{it:Note:} Gray areas denote 95 % confidence intervals", size(vsmall) span)) /// ytitle("% of {bf:women} reporting (suspected) infertility") xtitle("Age group") /// xlabel(1/5, val alternate) /// ysize(8) name(infertage, replace)
Oct 15, 2014
Random graphs (33): Regression parameters
use V4 V56 SEX using "ZA5900_v2-0-0.dta", replace // ISSP 2012 data
renvars, lower // Switch variable names to lower case
kountry v4, from(iso3n) to(iso2c) // Generate country abbreviation variable
replace _ISO2C_ = "RU" if v4 == 643
encode _ISO2C_, gen(country)
drop if cntry == "ZA" // v56 missing in ZA
// Prepare variables
recode sex (9 = .)
recode v56 (0 8 9 = .), gen(jobsat_rev)
generate jobsat = 7 - jobsat_rev
label define jobsat 6 "Completely satisfied" ///
3 "Neither satisfied nor dissatisfied" ///
0 "Completely dissatisfied"
label value jobsat jobsat
preserve
// Save parameters
statsby gendergap = _b[sex] ///
loci = (_b[sex] - 1.96 * _se[sex]) ///
hici = (_b[sex] + 1.96 * _se[sex]) ///
, by(country) clear total: ///
regress jobsat sex
replace country = 1000 if country == . // Label parameter from total sample
label define country 1000 "{bf: Total}", modify // Label parameter from total sample
egen order_ = rank(-gendergap), unique
labmask order_, value(country) decode
twoway (dot gendergap order_, msize(vsmall)) ///
(rcap loci hici order_) ///
, legend(off) xlabel(1/37, valuelabels ang(v)) ///
yline(0) note("{it:Source:} ISSP 2012 (Family and Changing Gender Roles IV), own calculations", span) ///
ylabel(,format(%6.1f)) name(gendergap, replace) ///
xtitle(" ") ytitle("Gender gap in job satisfaction")
restore
Oct 7, 2014
Random graphs (32): Dot plot with confidence intervals
// Read in data points
input str50 variable_str coeff se
"Discrimination in being hired" .18 .08
"Vocational index" 1.97 .53
end
// Switch to numeric
encode variable_str, gen(variable)
// Create variable for -by-
recode variable (1 = 2 "Model 2") (2 = 1 "Model 1"), gen(model)
// Calculate 95 % CI's
gen loci = coeff - 1.96 * se
gen hici = coeff + 1.96 * se
twoway (dot coeff variable, by(model, legend(off) note("")) ///
horizontal dotextend(no)) ///
(rcap hici loci variable, horizontal) ///
, legend(off) ylabel(1/2, valuelabels) ///
xline(0) ytitle("") ///
ysize(4) yscale(range(0 3))
Jun 20, 2014
Random graphs (26): Plot means by country
use cntry stflife using "ESS1e06.3_F1.dta", clear
encode cntry, gen(country) // De-string country identifier
preserve
statsby mean_ = _b[_cons] ///
loci = (_b[_cons] - 1.96 * _se[_cons]) ///
hici = (_b[_cons] + 1.96 * _se[_cons]) ///
, by(country) total clear: /// // -total- fits model for total sample
regress stflife
replace country = 1000 if country == . // Label mean from total sample
label define country 1000 "{bf: Total}", modify // Label mean from total sample
egen order_ = rank(-mean_), unique
labmask order_, value(country) decode
twoway (dot mean_ order_, msize(vsmall)) ///
(rcap loci hici order_) ///
, legend(off) xlabel(1/23, valuelabels ang(v)) ///
ylabel(,format(%6.1f)) name(lifesat, replace) ///
xtitle("") ytitle("Life satisfaction (avg.)")
restore
Jun 2, 2014
Random graphs (21): Confidence interval plots
clear
input str33 fiel prop_women prop_women_se ci_low ci_hi str3 test_
"Education" 0.7365 0.0029 0.7308 0.7422 No
"Social and behavioral sciences" 0.5251 0.0163 0.4931 0.5569 Yes
"Psychology" 0.8282 0.0183 0.7921 0.8639 No
"Nursing and caring" 0.9298 0.0046 0.9208 0.9388 No
"Therapy and rehabilitation" 0.8598 0.0336 0.7939 0.9257 No
"Child care and youth service" 0.9517 0.0178 0.9168 0.9866 No
"Social work and counseling" 0.8419 0.0238 0.7953 0.8885 No
end
// Transform from proportion to percentage
replace prop_women = prop_women * 100
replace ci_low = ci_low * 100
replace ci_hi = ci_hi * 100
// Sort educational fields by % females
egen order_ = rank(-prop_women), unique
labmask order_, val(fiel)
twoway (dot prop_women order_, horizontal) ///
(rcap ci_low ci_hi order_, horizontal) ///
, legend(off) ylabel(1/7, valuelabels) ///
ytitle("Educational fields") ///
xtitle("% Females") ///
xline(50) xlabel(40 (10) 100, format(%6.0f)) ///
note("{it: Source:} European Labor Force Survey 2009" ///
"{it: Note:} Error bars denote 95% CI's", span)
clear
input str44 fiel prop_women prop_women_se ci_low ci_hi hard
"Science, mathematics, and computing" 0.445 0.0114 0.4226 0.4673 1
"Physical science" 0.383 0.0069 0.3694 0.3964 1
"Physics" 0.390 0.0288 0.3336 0.4464 1
"Mathematics" 0.571 0.0193 0.5331 0.6089 1
"Computer science" 0.252 0.0051 0.1460 0.3580 1
"Engineering, manufacturing, and construction" 0.161 0.0012 0.1586 0.1634 1
"Engineering and engineering trades" 0.199 0.0182 0.1633 0.2347 1
"Humanities" 0.717 0.0096 0.698 0.7358 0
"Social and behavioral science" 0.525 0.0163 0.4931 0.5569 0
"Psychology" 0.828 0.0183 0.7921 0.8639 0
"Sociology and cultural studies" 0.663 0.0354 0.5936 0.7323 0
end
// Transform from proportion to percentage
replace prop_women = prop_women * 100
replace ci_low = ci_low * 100
replace ci_hi = ci_hi * 100
// Sort educational fields by % females and by hard/soft
egen orderhard = rank(-prop_women) if hard, unique
labmask orderhard, val(fiel)
egen ordersoft = rank(-prop_women) if hard == 0, unique
labmask ordersoft, val(fiel)
twoway (dot prop_women orderhard, horizontal) ///
(rcap ci_low ci_hi orderhard, horizontal) ///
, legend(off) ylabel(1/7, valuelabels) ///
ytitle("Hard educational fields") ///
xtitle("% Females") ///
xline(50) xlabel(20 (10) 100, format(%6.0f)) ///
xscale(off) /// // remove x-axis
name(hard, replace)
twoway (dot prop_women ordersoft, horizontal) ///
(rcap ci_low ci_hi ordersoft, horizontal) ///
, legend(off) ylabel(1/4, valuelabels) ///
ytitle("Soft educational fields") ///
xtitle("% Females") ///
graphregion(margin(l=28)) /// Account for different y-axis label length
xline(50) xlabel(20 (10) 100, format(%6.0f)) ///
name(soft, replace)
graph combine hard soft, col(1) imargin(b = 2 t = 1) ///
note("{it: Source:} European Labor Force Survey 2009" ///
"{it: Note:} Error bars denote 95% CI's", span)
graph export Graph.png
Mar 17, 2014
Random graphs (19): Simple interaction plots with -marginsplot- for continuous variables
sysuse auto, clear
regress mpg c.weight##c.gear_ratio
// Obtain means and standard deviation for plots
qui sum weight
local w_minsd = r(mean) - r(sd)
local w_mean = r(mean)
local w_plusd = r(mean) + r(sd)
qui sum gear_ratio
local g_minsd = r(mean) - r(sd)
local g_mean = r(mean)
local g_plusd = r(mean) + r(sd)
// Calculate predicted values for plot
margins, at(c.gear_ratio = (`g_minsd' `g_mean' `g_plusd') ///
c.weight = (`w_minsd' `w_mean' `w_plusd')) vsquish
// Plot
marginsplot, recastci(rarea) ciopts(color(gs10)) ///
title("Weight moderates the gear ratio{char 0150}mileage association") ///
xlabel(`g_minsd' "-1 SD" `g_mean' "Average gear ratio" `g_plusd' "+1 SD") ///
ytitle("Predicted mileage") ///
xtitle("Gear ratio") ///
plotopts(msymbol(none)) /// // Turn off markers
plot1opts(lpattern(longdash)) /// // Define line types here
plot2opts(lpattern(solid)) ///
plot3opts(lpattern(shortdash)) ///
legend(subtitle(Weight) ///
order(4 "- 1 SD" 5 "Mean" 6 "+ 1 SD" 2 "95% CI"))
// For some reason, marginsplot ignores the -label option-, therefore
// -order- is used here
gr export plot.png, as(png) replace
Feb 5, 2014
Random graphs (18): Simple interaction plots with -marginsplot- for categorical and continuous variables
sysuse auto, clear
reg mpg c.weight##i.foreign
qui sum weight
local minsd = r(mean) - r(sd)
local mean = r(mean)
local plusd = r(mean) + r(sd)
margins foreign, at(weight = (`minsd' `mean' `plusd'))
marginsplot, recast(line) noci scheme(lean1) ///
name(first, replace)
gr export first.png, as(png) replace
// Impossible to tell the difference between
// domestic and foreign lines because of the lean scheme
marginsplot, recast(line) noci ///
plot1opts(lpattern(dash)) ///
scheme(lean1) ///
name(second, replace)
gr export second.png, as(png) replace
// That fixes the problem
marginsplot, recast(line) recastci(rarea) ciopts(color(gs10)) ///
plot1opts(lpattern(dash)) ///
scheme(lean1) ///
name(third, replace)
gr export third.png, as(png) replace
// That's a nice way of plotting CI's
marginsplot, recast(line) recastci(rarea) ciopts(color(gs10)) ///
plot1opts(lpattern(dash)) ///
ytitle("Predicted mileage") ///
title("Car type moderates the weight{char 0150}mileage association") ///
xlabel(`minsd' "-1 SD" `mean' "Average weight (lbs.)" `plusd' "+1 SD") ///
xtitle("") ///
scheme(lean1) ///
name(fourth, replace)
gr export fourth.png, as(png) replace
// That labels the graph nicely
Oct 8, 2013
Random graphs (16): Plots with confidence intervals
preserve
reg srh i.country [pweight = weight]
predict means, xb
predict semeans, stdp
gen loci = means - 1.96 * semeans
gen hici = means + 1.96 * semeans
keep if pickone
egen order_ = rank(-means), unique
labmask order_, value(country) decode
twoway (dot means order_) (rcap loci hici order_) ///
, legend(off) xlab(1/29, valuelabels ang(v)) ///
ylabel(, format(%6.1f)) ///
ytitle("Self-rated health") ///
xtitle("") ///
note(`"Note: SRH ranges from 0 ("poor") to 4 ("excellent")"', span)
restore
preserve logit srh_bin i.country [pweight = weight] predict lr_index, xb predict se_index, stdp generate p_hat = exp(lr_index)/(1+exp(lr_index)) // Predicted probabilities // (following http://www.stata.com/support/faqs/...) gen lb = lr_index - invnormal(0.975)*se_index gen ub = lr_index + invnormal(0.975)*se_index gen plb = exp(lb)/(1+exp(lb)) gen pub = exp(ub)/(1+exp(ub)) replace p_hat = p_hat * 100 replace plb = plb * 100 replace pub = pub * 100 keep if pickone egen order_ = rank(p_hat), unique labmask order_, value(country) decode twoway (dot p_hat order_) (rcap plb pub order_) /// , legend(off) xlab(1/29, valuelabels ang(v)) /// ylabel(, format(%6.0f)) /// ytitle("% Poor health") xtitle("")
Mar 31, 2013
Random graphs (13): Interaction plots using -marginsplot-
estimates use m2 // Open estimates saved via -estimates save-
xtcloglog // Show results
// Set esample() as explained in the help-file to -estimates save-:
mark miss_mod1
markout miss_mod1 fail2 female age1519 age2529 age3034 ///
hatlev1d loggdp hatfield_main
estimates esample: if miss_mod1
// Calculate predicted probabilities:
margins female#hatlev1d, pred(pu0)
// Simple default marginsplot:
marginsplot, xtitle("Gender") ytitle("Predicted probability") ///
caption("{it:Source:} EU-LFS 2009 AHM", span) ///
title("") name(sexXeduc1, replace)
// As a bar graph:
marginsplot, xtitle("") ytitle("Predicted probability") title("") ///
recast(bar) /// // Draw bar graph instead of default
xdimension(female hatlev1d) ///
/* Get both interacting variables onto x-axis */ ///
caption("{it:Source:} EU-LFS 2009 AHM", span) ///
name(sexXeduc2, replace)
// As a bar graph that better separates the different variable combinations:
marginsplot, xtitle("") ytitle("Predicted probability") ///
recast(bar) ///
xdimension(female) ///
bydimension(hatlev1d, ///
elabels(1 "Low education" ///
2 "Medium education" ///
3 "High education")) ///
/* Use -bydimension- to better structure the x-axis */ ///
name(sexXeduc3, replace)
// As a proper-looking graph:
marginsplot, xtitle("") ytitle("Predicted probability") ///
recast(bar) ///
xdimension(female) ///
bydimension(hatlev1d, elabels(1 "Low education" ///
2 "Medium education" ///
3 "High education")) ///
byopts(row(1) /// All plots in one row
noiyaxes /// // Don't display individual y-axes
imargin(zero) /// // No margin between by() dimensions
title("") /// // Suppress title
caption("{it:Note:} Error bars are 95 % CI's" ///
"{it:Source:} EU-LFS 2009 AHM", span)) ///
/* Caption here to suppress displaying it three times */ ///
plotregion(lwidth(none)) ///
/* Removes lines around separate by dimension plots */
subtitle(, pos(6)) /// // Place label of by dimensions below plot
plotopts(fcolor(gs14) lcolor(b)) ///
/* Options of the recast() type can go here as well */ ///
name(sexXeduc4, replace)
Mar 18, 2013
Random graphs (8): Coefficient plot by hand
// Requires parmest, ingap, and eclplot
use eulfsahm2010, clear
// Create and label dummy variables for regression model
fre age
quietly tab age, gen(age)
label var age1 "15-19 years"
label var age2 "20-24 years"
label var age3 "25-29 years"
label var age4 "30-34 years"
label var age5 "35-39 years"
label var age6 "40-44 years"
label var age7 "45-49 years"
label var age8 "50-54 years"
label var age9 "55-59 years"
label var age10 "60-64 years"
fre marstat
quietly tab marstat, gen(marstat)
label var marstat1 "Widowed, divorced, or separated"
label var marstat2 "Single"
label var marstat3 "Married"
fre education
quietly tab education, gen(education)
label var education1 "Low education"
label var education2 "Medium education"
label var education3 "Tertiary education"
fre parttime
quietly tab parttime, gen(parttime)
label var parttime1 "Marginal part-time (< 20 hrs.)"
label var parttime2 "Substantial part-time (20-34 hrs.)"
label var parttime3 "Full-time (35-80 hrs.)"
fre occgroups
quietly tab occgroups, gen(occgroups)
label var occgroups1 "Unskilled job"
label var occgroups2 "Skilled manual job"
label var occgroups3 "Skilled routine services job"
label var occgroups4 "High-skilled services"
fre industry
quietly tab industry , gen(industry)
label var industry1 "Agriculture"
label var industry2 "Industry and construction"
label var industry3 "Market services"
label var industry4 "Non-market services"
// Re-label some more variables here.
// Variable labels HAVE to be identical to those specified below.
// Formatting is possible here.
label var female "{bf:Female sex}"
label var smallfirm "{bf:Small firm} (< 10 employees)"
label var fixedterm "{bf:Fixed-term contract}"
label var jobtenure "{bf:Job tenure} (in years)"
label var supervisor "{bf:Supervisor status}"
// Estimate regression model
ologit posstend female age1-age3 age5-age10 marstat1 marstat3 ///
education1 education3 parttime1 parttime2 smallfirm fixedterm ///
jobtenure supervisor occgroups1 occgroups3 occgroups4 ///
industry1 industry3 industry4, cluster(cntry)
// Save regression parameters in data set
parmest, norestore eform label
// Add gaps to add some structure to the regression plots
ingap 2 11 13 15 21 24, rowlabel(label) ///
growlabel("Age ({it:Ref.} 30-34 years)" ///
"Marital status ({it:Ref.} Single)" ///
"Education ({it:Ref.} Medium)" ///
"Working hours ({it:Ref.} Full-time (35-80 hrs)" ///
"Occup. group ({it:Ref.} Skilled manual job)" ///
"Industry ({it:Ref.} Industry and construction)")
// Define label by hand to be able to determine order of predictors on axis
// Labels must be exactly the same as above
//
// Formatting axis label text seems
// to be working here as well
//
// An easier solution might be the one here:
// http://www.survey-design.com.au/Stata%20Graphs.html
// ("Odds ratio graph - selecting order of bars")
#delimit ;
label define label1 1 "{bf:Female sex}"
2 "{bf:Age} ({it:Ref.} 30-34 years)"
3 "15-19 years"
4 "20-24 years"
5 "25-29 years"
6 "35-39 years"
7 "40-44 years"
8 "45-49 years"
9 "50-54 years"
10 "55-59 years"
11 "60-64 years"
12 "{bf:Marital status} ({it:Ref.} Single)"
13 "Widowed, divorced, or separated"
14 "Married"
15 "{bf:Education} ({it:Ref.} Medium)"
16 "Low education"
17 "Tertiary education"
18 "{bf:Working hours} ({it:Ref.} Full-time (35-80 hrs)"
19 "Marginal part-time (< 20 hrs.)"
20 "Substantial part-time (20-34 hrs.)"
21 "{bf:Small firm} (< 10 employees)"
22 "{bf:Fixed-term contract}"
23 "{bf:Job tenure} (in years)"
24 "{bf:Supervisor status}"
25 "{bf:Occupational group} ({it:Ref.} Skilled manual job)"
26 "Unskilled job"
27 "Skilled routine services job"
28 "High-skilled services"
29 "{bf:Industry} ({it:Ref.} Industry and construction)"
30 "Agriculture"
31 "Market services"
32 "Non-market services"
33 "Constant"
34 "Constant" ;
#delimit cr
encode label, gen(label1) label(label1)
// Specifying a label in the -encode- command makes sure that -encode-
// does not sort labels in alphabetical order
eclplot estimate min95 max95 label1 ///
if label1 <= 32 /// Don't plot the cut values
, horizontal ylabel(1(1)32, labsize(vsmall)) /// Reduce font size on axis
xline(1) /// Add vertical line
xtitle("Odds ratio") ///
ytitle("Predictors") ///
title("Possibility to vary start and/or end of working day", span) ///
caption("Source: EU-LFS AHM 2010, authors' calculations." /*
*/ "Ordered logit regression, /*
*/ 95% CI's adjusted for clustering in countries", span)
Nov 16, 2012
Random Graphs (5): Combining plots of means
set autotabgraphs on // Recycle graph window
import excel "X:\Data\Desktop\means.xlsx", sheet("Sheet1") ///
firstrow clear case(lower)
// Read in aggregate data from Excel sheet:
// - Variable name (variable) [string]
// - Country name (country) [string]
// - Average score (mean)
// - Confidence bounds (lerr uerr)
// Create numerical country variable for x-axis:
// Fast alternative: Use -encode- if you do not
// want to manipulate the order of countries
generate cntry = 1 if country == "SE"
replace cntry = 2 if country == "UK"
replace cntry = 3 if country == "NL"
replace cntry = 4 if country == "DE"
replace cntry = 5 if country == "PT"
label define cntry 1 "SE" 2 "UK" 3 "NL" 4 "DE" 5 "PT"
label val cntry cntry
label var cntry "Country"
graph drop _all // Drop graphs in memory, usefulfor re-running this file
twoway (scatter mean cntry if variable == "Job autonomy") ///
(rcap lerr uerr cntry if variable == "Job autonomy") ///
, legend(off) xlabel(, valuelabels) ///
ytitle("Job autonomy") xscale(off) name(JA)
twoway (scatter mean cntry if variable == "WF culture") ///
(rcap lerr uerr cntry if variable == "WF culture") ///
, legend(off) xlabel(, valuelabels) ///
ytitle("WF culture") xscale(off) name(WFC)
twoway (scatter mean cntry if variable == "WF supervisor support") ///
(rcap lerr uerr cntry if variable == "WF supervisor support") ///
, legend(off) xlabel(, valuelabels) ///
ytitle("WF supervisor support") xscale(off) name(WFSS)
twoway (scatter mean cntry if variable == "WF co-worker support") ///
(rcap lerr uerr cntry if variable == "WF co-worker support") ///
, legend(off) xlabel(, valuelabels) ///
ytitle("WF co-worker support") xscale(off) name(WFCS)
twoway (scatter mean cntry if variable == "FWA use") ///
(rcap lerr uerr cntry if variable == "FWA use") ///
, legend(off) xlabel(, valuelabels) ///
ytitle("FWA use") name(FWA)
graph combine JA WFC WFSS WFCS FWA, col(1) xcommon imargin(b = 0 t = 0) ysize(10)
Subscribe to:
Posts (Atom)