clear
input str100 grouping t1 t2
"Department managers" .22 .399
"Less than year" .35 .21
"Remaining non-managers" .31 .36
end
reshape long t, i(grouping) j(treatment)
twoway (connect t treatment if grouping == "Department managers", msymbol(o)) ///
(connect t treatment if grouping == "Less than year", msymbol(o)) ///
(connect t treatment if grouping == "Remaining non-managers", msymbol(o)) ///
(scatteri .399 2 "Department managers" ///
.21 2 "Less than year" ///
.36 2 "Remaining non-managers", msymbol(none)) ///
, legend(off) xtitle("Minimum wage raise") xlabel(0 " " 1 "Before" 2 "After" 3 " ", notick) ///
ytitle("Job satisfaction") ylabel(, format(%6.2f))
Showing posts with label reshape. Show all posts
Showing posts with label reshape. Show all posts
Jan 7, 2018
Random graphs (124): Line plots
Labels:
Random graphs,
reshape,
twoway connected,
twoway scatteri
Nov 17, 2017
Random graphs (119): Line plots and dot plots
// Prepare ESS round 8
use cntry essround wrkctra dweight using "ESS8e01.dta", clear
recode wrkctra (6 = .a) (7 = .b) (8 = .c) (9 = .d)
// Add ESS rounds 1-7
append using "ESS1-7e01.dta", keep(cntry essround wrkctr wrkctra dweight)
// Prepare variables
generate nocontract = (wrkctra == 3) if !missing(wrkctra)
kountry cntry, from(iso2c)
rename NAMES_STD country
// Prepare files for post commands
tempname foo1
tempname foo2
postfile `foo1' str20 country essround ll nocontract ul using `foo2', replace
qui levelsof country, local(country)
// Loop
foreach x of local country {
foreach i of numlist 1/8 {
capture logit nocontract [pw = dweight] if country == "`x'" & essround == `i'
if _rc == 0 {
qui margins [pw = dweight]
matrix prevs = r(table)
local prev = prevs[1,1] * 100
local ll = prevs[5,1] * 100
local ul = prevs[6,1] * 100
*di "`x'" _skip(5) `i' _skip(5) `ll' _skip(5) `prev' _skip(5) `ul'
post `foo1' ("`x'") (`i') (`ll') (`prev') (`ul')
}
}
}
postclose `foo1'
// Plot estimates
use `foo2', clear
// Fix value label
label define essround 1 "2002" 2 "2004" 3 "2006" 4 "2008" ///
5 "2010" 6 "2012" 7 "2014" 8 "2016", modify
label val essround essround
// First plot
sort country essround
twoway (rarea ll ul essround, lcolor(white)) ///
(connected nocontract essround), ///
by(country, ///
note("") ///
legend(off)) ///
xlabel(1/8, val ang(90)) ///
xtitle("") ytitle("% without job contract") ///
name(figure2a, replace) ysize(9)
// Second plot
reshape wide ll ul nocontract, i(country) j(essround)
scores average = mean(nocontract*)
egen order_ = rank(average), unique
labmask order_, value(country)
twoway (dot nocontract1 order_, horizontal) ///
(dot nocontract2 order_, horizontal) ///
(dot nocontract3 order_, horizontal) ///
(dot nocontract4 order_, horizontal) ///
(dot nocontract5 order_, horizontal) ///
(dot nocontract6 order_, horizontal) ///
(dot nocontract7 order_, horizontal) ///
(dot nocontract8 order_, horizontal) ///
(rspike ul1 ll1 order_, horizontal) ///
(rspike ul2 ll2 order_, horizontal) ///
(rspike ul3 ll3 order_, horizontal) ///
(rspike ul4 ll4 order_, horizontal) ///
(rspike ul5 ll5 order_, horizontal) ///
(rspike ul6 ll6 order_, horizontal) ///
(rspike ul7 ll7 order_, horizontal) ///
(rspike ul8 ll8 order_, horizontal), ///
ylabel(1/32, val) ytitle("") xscale(alt) ///
xtitle("% without job contract") ///
legend(order(1 "2002" 2 "2004" 3 "2006" 4 "2008" ///
5 "2010" 6 "2012" 7 "2014" 8 "2016") ///
pos(5) ring(0)) ///
name(figure2b, replace) ysize(9)
graph combine figure2a figure2b, ///
col(2) name(figure2, replace) ///
note(" " "{it:Source:} European Social Survey 2002-16, weighted data. {it:Note:} Error bands/spikes denote 95% confidence intervals", span size(*.7))
Labels:
European Social Survey,
foreach,
graph combine,
kountry,
labmask,
logit,
margins,
postfile,
Random graphs,
reshape,
scores,
twoway dot,
twoway rarea,
twoway rspike,
weight
Jul 27, 2017
Random graphs (105): Scatterplot
// Open ISSP data
use C_ALPHAN SEX V51 WEIGHT using ZA4350_v2-0-0.dta, clear
// Fix country variable
replace C_ALPHAN = "BE" if C_ALPHAN == "BE-FLA"
replace C_ALPHAN = "DE" if C_ALPHAN == "DE-E"
replace C_ALPHAN = "DE" if C_ALPHAN == "DE-W"
replace C_ALPHAN = "UK" if C_ALPHAN == "GB-GBN"
// Recode job satisfaction
recode V51 (7 = 0 "Completely dissatisfied") ///
(6 = 1 "Very dissatisfied") ///
(5 = 2 "Fairly dissatisfied") ///
(4 = 3 "Neither satisfied nor dissatisfied") ///
(3 = 4 "Fairly satisfied") ///
(2 = 5 "Very satisfied") ///
(1 = 6 "Completely satisfied"), gen(jobsat)
label var jobsat "Job satisfaction"
// Recode sex
gen female = (SEX == 2) if !missing(SEX)
label define female 0 "Male" 1 "Female"
label val female female
label var female "Sex"
drop if missing(female)
// Get rid of old variables
drop SEX V51
// Collapse and reshape data
collapse jobsat [pw = WEIGHT], by(C_ALPHAN female)
reshape wide jobsat, i(C_ALPHAN) j(female)
label var jobsat0 "Job satisfaction: men"
label var jobsat1 "Job satisfaction: women"
// Plot
twoway (scatter jobsat0 jobsat1 if jobsat0 > jobsat1, mlabel(C_ALPHAN) mlabpos(9)) ///
(scatter jobsat0 jobsat1 if jobsat0 < jobsat1, mlabel(C_ALPHAN) mlabpos(6)) ///
(function y = x, range(3.5 5.1)) ///
, xlabel(3.5 (.5) 5.0, format(%6.2f)) ///
ylabel(3.5 (.5) 5.0, format(%6.2f)) ///
ytitle("Job satisfaction: men") ///
xtitle("Job satisfaction: women") ///
title("Cross-national variation in the gender gap in job satisfaction", span) ///
note(" " "{it:Source:} International Social Survey Program (ISSP) 2005, doi:10.4232/1.11648", span) ///
legend(off) name(figure1, replace)
Feb 26, 2017
Chapter 2 of Singer and Willett's (2003) book on longitudinal data analysis
There is a different take on the chapter here, but I like mine better.
// Open data set use http://www.ats.ucla.edu/stat/stata/examples/alda/data/tolerance, clear // Figure 2.1 list, sep(0) noobs reshape long tol, i(id) j(age) list if inlist(id, 9, 45, 1653), sepby(id) noobs // Table 2.1 reshape wide quietly estpost cor tol*, matrix listwise esttab, unstack not noobs compress nostar b(%6.2f) // Figure 2.2
reshape long twoway (scatter tol age), by(id) ytitle(Tolerance) xtitle(Age) name(Figure22, replace) // Figure 2.3
twoway (scatter tol age) ///
(lowess tol age) ///
, by(id, note("Graphs by id, lowess curves") legend(off)) ///
ytitle(Tolerance) xtitle(Age) name(Figure23, replace)
// Table 2.2
preserve
tempfile table22
generate time = age - 11
statsby initial = _b[_cons] ///
initial_se = _se[_cons] ///
change = _b[time] ///
change_se = _se[time] ///
residual = (e(rmse)^2) ///
explained = e(r2) ///
, by(id) saving(`table22', replace): regress tol time
drop time
reshape wide
merge 1:1 id using `table22', nogenerate
format initial initial_se change change_se residual explained exposure %6.2f
keep id initial initial_se change change_se residual explained male exposure
order id initial initial_se change change_se residual explained male exposure
list, sep(0) noobs abbrev(15)
tempfile figure28
save `figure28', replace
// Figure 2.4
stem initial, round(.01)
stem change, round(.01)
stem residual, round(.01)
stem explained, round(.01)
restore
// Figure 2.5
twoway (scatter tol age) ///
(lfit tol age) ///
, by(id, note("Graphs by id, OLS curves") legend(off)) ///
ytitle(Tolerance) xtitle(Age) name(Figure25, replace)
// Figure 2.6
xtset id age
xtline tol, overlay t(age) i(id) legend(off) ///
ytitle(Tolerance) ylabel(1 (1) 4) ///
addplot(lowess tol age, lwidth(thick) lpattern(solid)) ///
xtitle(Age) title("Observed data and lowess smoother") ///
xsize(3) name(Figure26A, replace) nodraw
quietly regress tol i.id##c.age
predict tolhat
xtline tolhat, overlay t(age) i(id) legend(off) ///
ytitle(Tolerance) ylabel(1 (1) 4) ///
addplot(lfit tol age, lwidth(thick) lpattern(solid)) ///
xtitle(Age) title("OLS trajectories") ///
xsize(3) name(Figure26B, replace) nodraw
graph combine Figure26A Figure26B, col(2) name(Figure26, replace)
// Table 2.2
preserve
use `table22', clear
tabstat initial change, stat(mean sd) format(%6.2f)
cor initial change
restore
// Figure 2.7
quietly sum exposure, detail
generate highexposure = (exposure >= r(p50))
generate lowexposure = (exposure < r(p50))
generate female = !male
label var male "Males"
label var female "Females"
label var lowexposure "Low exposure"
label var highexposure "High exposure"
foreach z of varlist female male lowexposure highexposure {
capture drop tolhat
quietly regress tol i.id##c.age if `z'
predict tolhat if `z'
xtline tolhat if `z', ///
overlay t(age) i(id) legend(off) ///
ytitle("Predicted tolerance") ylabel(1 (1) 4) ///
addplot(lfit tol age if `z', lwidth(thick) lpattern(solid)) ///
xtitle(Age) title(`: variable label `z'') ///
xsize(3) ysize(2) name(Figure27_`z', replace) nodraw
}
graph combine Figure27_female Figure27_male Figure27_lowexposure Figure27_highexposure, col(2) name(Figure27, replace)
// Figure 2.8
use `figure28', clear
quietly correlate initial male
local pmca = round(r(rho), .01)
quietly correlate initial exposure
local pmcb = round(r(rho), .01)
quietly correlate change male
local pmcc = round(r(rho), .01)
quietly correlate change exposure
local pmcd = round(r(rho), .01)
label define male -1 " " 1 "Male" 0 "Female" 2 " ", modify
label val male male
twoway dot initial male, ///
ytitle("Predicted intercept") ///
xtitle("Gender") ///
xlabel(-1 0 1 2, val) ///
ysize(2) xsize(3) ///
note("{it:r} = `pmca'") ///
name(Figure28A, replace) nodraw
twoway scatter initial exposure, ///
ytitle("Predicted intercept") ///
xtitle("Exposure") ///
ysize(2) xsize(3) ///
note("{it:r} = `pmcb'") ///
name(Figure28B, replace) nodraw
twoway dot change male, ///
ytitle("Predicted change") ///
xtitle("Gender") ///
ysize(2) xsize(3) ///
xlabel(-1 0 1 2, val) ///
note("{it:r} = `pmcc'") ///
name(Figure28C, replace) nodraw
twoway scatter change exposure, ///
ytitle("Predicted change") ///
xtitle("Exposure") ///
ysize(2) xsize(3) ///
note("{it:r} = `pmcd'") ///
name(Figure28D, replace) nodraw
graph combine Figure28A Figure28B Figure28C Figure28D, col(2) row(2) name(Figure28, 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
Labels:
estpost,
reshape,
statsby,
Textbooks,
twoway dot,
twoway scatter,
xtline,
xtset
Apr 25, 2016
Random graphs (74): Arrow plot
use 7363_F1.dta, clear
// Prepare wave variable
label define wavex 1 "1991" 2 "1995" 3 "2000" 4 "2005" 5 "2010"
label val wave wavex
// Using the internet
fre y10_q_16
tab y10_q_16 wave
keep if wave > 3 // Drop waves where variable is missing
// Create variable of interest
generate nointernetuse = .
replace nointernetuse = (y10_q_16 == 7) if !missing(y10_q_16)
// Collapse data
collapse (mean) nointernetuse, by(countid wave)
replace nointernetuse = nointernetuse * 100
reshape wide nointernetuse, i(countid) j(wave)
// Sort variable labels
egen order_ = rank(-nointernetuse5), unique
labmask order_, val(countid) decode
label define order_ 8 "Macedonia", modify
twoway (pcarrow nointernetuse4 order_ nointernetuse5 order_, horizontal) ///
(scatter order_ nointernetuse5 if missing(nointernetuse4)) ///
, ylabel(1/34, val) legend(order(1 "Decrease 2005 to 2010" ///
2 "No data for 2005") pos(7) ring(0) size(*.8)) ///
ytitle("") ysize(8) xscale(alt) xtitle("% of labor force" ///
"never using internet/email" ///
"for professional reasons") ///
xlabel(10 (10) 80, grid) ///
note(" " "{it:Source:} European Working Conditions Survey 2005 and 2010", span)
Apr 13, 2016
Random graphs (72): Dot plots with the by() option
clear
input str20 cntry bindirect blbindirect ulbindirect bdirect blbdirect ulbdirect
Austria 0.002 -0.002 0.005 0.008 -0.016 0.032
Germany 0.006 0.000 0.012 0.040 0.017 0.063
Sweden 0.001 -0.005 0.007 0.054 0.031 0.076
Netherlands 0.000 -0.006 0.006 0.039 0.017 0.061
Spain 0.005 0.002 0.009 0.016 -0.001 0.033
Italy 0.000 -0.005 0.005 0.029 0.011 0.047
France 0.001 -0.003 0.006 0.043 0.025 0.061
Denmark 0.001 0.000 0.003 0.015 0.001 0.029
Greece 0.000 -0.003 0.003 0.041 0.026 0.057
Switzerland 0.002 -0.002 0.006 0.003 -0.020 0.026
Belgium 0.000 -0.003 0.004 0.033 0.016 0.050
"Czech Republic" 0.004 -0.001 0.010 0.073 0.031 0.115
Poland -0.004 -0.010 0.002 0.034 -0.003 0.071
end
// Reshape to long format
reshape long b blb ulb, j(fx "indirect" "direct") i(cntry) string
// Encode variables
encode cntry, gen(country)
encode fx, gen(effect)
label define effect 1 "Direct effect" 2 "Indirect effect", modify
// Sort by biggest direct effect size
egen order_ = rank(b) if effect == 1, unique
labmask order_, val(country) decode
bysort country (order_): replace order_ = order_[1] // Copy value to all cases
// xrescale makes sure that x-axis fits both by() plots
twoway (dot b order_, horizontal by(effect, legend(off) note(" ") xrescale)) ///
(rspike blb ulb order_, horizontal by(effect)) ///
, ylabel(1 (1) 13, val) xscale(alt) ///
ytitle(" ") xline(0)
Apr 11, 2016
Snijders and Bosker's (2012) chapter on multivariate multilevel analysis using -mixed-
Download the complete do-file including data here.
version 13.1
clear
// Data set for Example 16.1
// mlbook2, data set with pupils having missings on ses or IQ_verb or (langPOST and aritPOST) excluded.
input schoolnr pupilNR_new langPOST aritPOST ses IQ_verb IQ_perf Minority denomina sch_ses sch_iqv sch_min
1.00 1.00 . 12.00 -17.73 -1.37 -3.75 1.00 1.00 -14.035 -1.4039 0.630
// .. Skip a couple of thousand observations
end
// Rename outcome variables for -reshape-:
rename langPOST score1
rename aritPOST score2
// Convert data into long format:
reshape long score, i(pupilNR_new) j(outcome)
// Create dummy variables that identify outcome variable:
quietly tab outcome, gen(outcome)
// Results in Table 16.1:
mixed score outcome1 outcome2, nocons ///
|| schoolnr: outcome1 outcome2, nocons cov(un) ///
|| pupilNR_new: , nocons cov(un) residuals(un, t(outcome))
// Results in Table 16.2:
mixed score outcome1 outcome2 ///
outcome#c.IQ_verb outcome#c.ses ///
outcome#c.sch_iqv outcome#c.sch_ses ///
outcome#c.IQ_verb#c.ses outcome#c.sch_iqv#c.sch_ses, nocons ///
|| schoolnr: outcome1 outcome2, nocons cov(un) ///
|| pupilNR_new:, nocons cov(un) residuals(un, t(outcome))
Reference
Snijders, Tom, and Roel Boskers. 2012. Multilevel Analysis. An Introduction to Basic and Advanced Multilevel Modeling, 2nd ed. Sage.
Jan 5, 2016
Random graphs (55): Box plots
wbopendata, indicator(NY.GDP.PCAP.KD) clear
// Select countries
#delimit ;
keep if inlist(countrycode, "AUS", "AUT", "BEL", "CAN", "DNK") |
inlist(countrycode, "FIN", "FRA", "DEU", "GRC", "ISL") |
inlist(countrycode, "IRL", "ITA", "JPN", "NLD") |
inlist(countrycode, "NZL", "NOR", "PRT", "ESP") |
inlist(countrycode, "SWE", "CHE", "USA") ;
#delimit cr
// Get rid of unnecessary years
keep countryname countrycode iso2code region regioncode ///
indicatorname indicatorcode ///
yr1960 yr1970 yr1980 yr1990 yr2000 yr2010
// Reshape into long format
reshape long yr, i(countryname) j(year)
// Prepare outcome variable
rename yr gdppc
label var gdppc "GDP per capita (constant 2005 US$)"
// Plot
graph box gdppc, over(year)
Labels:
graph box,
inlist(),
Random graphs,
reshape,
wbopendata
Oct 4, 2015
Firebaugh (1997): Analyzing Repeated Surveys in Stata (ch. 4)
Note how this is restricted to the linear decompositions!
Intracohort change versus overall change
use year race age cohort racmar using GSS7214_R4.DTA, clear
// Select sample
keep if inlist(year, 1974, 1984, 1994)
keep if race == 1 // Whites only
keep if age >= 20
keep if !missing(age, racmar)
// Generate 10-year cohort variable
generate cohort10 = .
local lowerbracket = 1965
local upperbracket = 1974
foreach x of numlist 1/10 {
replace cohort10 = `x' if cohort >= `lowerbracket' & cohort <= `upperbracket'
label define cohort10 `x' "`lowerbracket' to `upperbracket'", modify
local lowerbracket = `lowerbracket' - 10
local upperbracket = `upperbracket' - 10
}
label value cohort10 cohort10
label var year "Year"
// Opposing interracial marriage
recode racmar (1 = 1 "Yes") (2 = 0 "No") (.d .i .n = .), gen(oppose)
label var oppose "Opposes interracial marriage"
generate oppose_perc = oppose * 100
// Create Table 4.1 without change
table cohort10 year, contents(mean oppose_perc n oppose_perc) format(%9.1f) center
// Different version of Table 4.1, including calculated changed
preserve
qui table cohort10 year, contents(mean oppose_perc n oppose_perc) format(%9.1f) center replace
rename table1 oppose_perc
rename table2 oppose_n
reshape wide oppose_perc oppose_n, i(cohort10) j(year)
gen change1984 = oppose_perc1984 - oppose_perc1974
gen change1994 = oppose_perc1994 - oppose_perc1984
format change1984 change1994 %9.1f
list, noobs sep(0) abbreviate(20) subvar
restore
Empirical example for linear decomposition: Trend in antiblack prejudice
use year race cohort racdin racpush racseg racmar using GSS7214_R4.DTA, clear
// Create sample
keep if inlist(year, 1972, 1976, 1980, 1984)
keep if race == 1
// Create outcome
alpha racdin racpush racseg racmar, item
egen zracdin = std(racdin)
egen zracpush = std(racpush)
egen zracseg = std(racseg)
egen zracmar = std(racmar)
scores y = total(zracdin zracpush zracseg zracmar), minvalid(3)
gen prejudice = 6 - y // Difficult to figure out the exact outcome
drop y // But this seems to be close enough
// Calculate total change in prejudice
mean prejudice, over(year) coeflegend
local totalchange = _b[1984] - _b[1972]
// Calculate annual change within cohorts
reg prejudice year cohort
local b1 = _b[year]
local b2 = _b[cohort]
// Calculate contributions of intracohort change and cohort replacement
local intracohort = `b1' * (1984 -1972)
di `intracohort'
qui mean cohort, over(year) coeflegend
local cohortreplacement = `b2' * (_b[1984] - _b[1972])
di `cohortreplacement'
// Calculate estimated total change
local estim_change = `intracohort' + `cohortreplacement'
di `estim_change'
// Output results
matrix results = (`totalchange')\(`b1')\(`b2')\(`intracohort')\(`cohortreplacement')\(`estim_change')
esttab matrix(results, fmt(%9.2f)), coeflabel(r1 "Total change" ///
r2 "Intracohort slope" ///
r3 "Intercohort slope" ///
r4 "Estimated contribution of intracohort change" ///
r5 "Estimated contribution of cohort replacement" ///
r6 "Estimated total change") ///
mtitle("Result") ///
varwidth(45)
Empirical example of the same-sign rule: Gender role attitudes
use year cohort fework fepres fepol fehome if year >= 1972 & year <= 1988 using GSS7214_R4.DTA, clear // Prepare variables recode fework fepres (2 = 0) recode fepol fehome (2 = 1) (1 = 0) // Create sample for Table 4.2 preserve keep if inlist(year, 1972, 1974, 1988) fre year gen id = _n reshape wide fework fepres fepol fehome, i(id) j(year) replace fepol1972 = fepol1974 replace fehome1972 = fehome1974 reshape long drop if year == 1974 recode year (1972 = 0 "1972") (1988 = 1 "1988"), gen(yr) // Table 4.2
foreach y of varlist fework fepres fepol fehome {
qui regress `y' yr
local mean1972 = _b[_cons]
local mean1988 = _b[_cons] + _b[yr]
local differen = _b[yr]
local t = _b[yr] / _se[yr]
*di "`y'" _skip(5) `mean1972' _skip(5) `mean1988' _skip(5) `differen' _skip(5) `t'
matrix stuff1 = (`mean1972' , `mean1988' , `differen' , `t')
if "`y'" == "fework" matrix table42 = stuff1
else matrix table42 = (table42\stuff1)
}
matrix rownames table42 = "WORK" "PRES" "POLI" "HOME"
matrix colnames table42 = "Mean 1972" "Mean 1988" "Change" "t-value"
esttab matrix(table42, fmt(3 3 3 1))
restore
// Table 4.3
foreach y of varlist fework fepres fepol fehome {
qui logit `y' year cohort
local n = e(N)
local intracohort = _b[year]
local intracohort_t = _b[year] / _se[year]
local intercohort = _b[cohort]
local intercohort_t = _b[cohort] / _se[cohort]
*di "`y'" _skip(5) `n' _skip(5) `mean1988' _skip(5) `differen' _skip(5) `t'
matrix stuff = (`n', `intracohort', `intracohort_t', ///
`intercohort', `intercohort_t')
if "`y'" == "fework" matrix table43 = stuff
else matrix table43 = (table43\stuff)
}
matrix rownames table43 = "WORK" "PRES" "POLI" "HOME"
matrix colnames table43 = "N" "Within cohort" "t" "Cross-cohort" "t"
esttab matrix(table43, fmt(0 %9.3f %9.1f %9.3f %9.1f)), modelwidth(14)
Reference
Firebaugh, Glenn. 1997. Analyzing Repeated Surveys. Sage. doi: 10.4135/9781412983396
Labels:
esttab,
foreach,
General Social Survey,
matrix,
matrix colname,
matrix rowname,
reshape,
table,
Textbooks
Aug 27, 2015
Preacher et al. (2008): Latent Growth Curve Modeling in Stata
The website accompanying Preacher et al.'s (2008) book on latent growth curve modeling provide syntax files for Lisrel, Mplus, and Mx, but not for Stata. These are the Stata commands for replicating the first four models presented in the book.
version 13 // Chapter 2 clear // Read in data from clsn_cov.dat ssd init clsn1 clsn3 clsn4 clsn5 clsn6 sex ssd set observations 851 ssd set means 37.9542 37.2785 37.0463 36.5696 36.1363 0.49 #delimit ; ssd set cov 6.3944 \ 3.2716 7.5282 \ 4.1435 6.0804 10.7290 \ 3.7058 5.1597 6.5672 10.2920 \ 4.1286 5.7608 7.2365 7.6463 12.9085 \ -0.0940 -0.0390 -0.1521 -0.1104 -0.1469 0.2502 ; #delimit cr // Table 2.2 ssd list // Model 0: The null model // Only mean intercept and residual variance are estimated
sem (Intercept@1 -> clsn1, ) /// // Set loadings to 1
(Intercept@1 -> clsn3, ) /// // for estimating constant
(Intercept@1 -> clsn4, ) ///
(Intercept@1 -> clsn5, ) ///
(Intercept@1 -> clsn6, ) ///
(clsn1 <- _cons@a, ) /// // Constrain means to be
(clsn3 <- _cons@a, ) /// // the same
(clsn4 <- _cons@a, ) ///
(clsn5 <- _cons@a, ) ///
(clsn6 <- _cons@a, ) ///
, latent(Intercept ) ///
cov(Intercept@0 /// Constrain intercept variance to 0
e.clsn1@b e.clsn3@b /// // Constrain residual variances
e.clsn4@b e.clsn5@b /// // to be the same
e.clsn6@b) ///
nocapslatent
estimates store m0
estat gof, stats(chi2 rmsea indices residuals)
// Non-normed fit index (NNFI) is called
// Tucker-Lewis index (TLI) in Stata
// Model 1: Random intercept model (Table 2.3)
// Only mean intercept, intercept variance, and residual variance are
// estimated
sem (Intercept@1 -> clsn1, ) /// // Set loadings to 1
(Intercept@1 -> clsn3, ) /// // for estimating constant
(Intercept@1 -> clsn4, ) ///
(Intercept@1 -> clsn5, ) ///
(Intercept@1 -> clsn6, ) ///
(clsn1 <- _cons@a, ) /// // Constrain means to be
(clsn3 <- _cons@a, ) /// // the same
(clsn4 <- _cons@a, ) ///
(clsn5 <- _cons@a, ) ///
(clsn6 <- _cons@a, ) ///
, latent(Intercept ) ///
cov(e.clsn1@b e.clsn3@b /// // Constrain residual variances
e.clsn4@b e.clsn5@b /// // to be the same
e.clsn6@b) ///
nocapslatent
estimates store m1
di 5.27 / (5.27 + 4.68) // Random intercept model allows calculating an ICC
estat gof, stats(chi2 rmsea indices residuals)
// Poor model fit according to all tests
// Likelihood ratio test:
lrtest m1 m0 // Massive improvement in fit for Model 1, though
// Model 2: Fixed intercept, fixed slope model (Table 2.4)
// Only mean intercept, intercept variance, and residual variance are
// estimated
sem (Intercept@1 -> clsn1, ) /// // Constrain paths to be 1
(Intercept@1 -> clsn3, ) ///
(Intercept@1 -> clsn4, ) ///
(Intercept@1 -> clsn5, ) ///
(Intercept@1 -> clsn6, ) ///
(clsn1 <- _cons@a, ) /// // Constrain intercepts to be the same
(clsn3 <- _cons@a, ) ///
(clsn4 <- _cons@a, ) ///
(clsn5 <- _cons@a, ) ///
(clsn6 <- _cons@a, ) ///
(Slope@0 -> clsn1, ) /// // Determine temporal structure
(Slope@2 -> clsn3, ) ///
(Slope@3 -> clsn4, ) ///
(Slope@4 -> clsn5, ) ///
(Slope@5 -> clsn6, ) ///
, covstruct(_lexogenous, diagonal) ///
latent(Intercept Slope ) ///
cov(Intercept@0 /// // Set intercept variance to 0
Slope@0 /// // Set slope variance to 0
e.clsn1@b ///
e.clsn3@b ///
e.clsn4@b ///
e.clsn5@b ///
e.clsn6@b) ///
means(Slope) /// // Estimate slope
nocapslatent
estimates store m2
estat gof, stats(chi2 rmsea indices residuals)
// Poor fit, even worse than Model 1
// Model 3: Random intercept, fixed slope (Table 2.5)
sem (Intercept@1 -> clsn1, ) /// // Constrain paths to be 1
(Intercept@1 -> clsn3, ) ///
(Intercept@1 -> clsn4, ) ///
(Intercept@1 -> clsn5, ) ///
(Intercept@1 -> clsn6, ) ///
(clsn1 <- _cons@a, ) /// // Constrain intercepts to be the same
(clsn3 <- _cons@a, ) ///
(clsn4 <- _cons@a, ) ///
(clsn5 <- _cons@a, ) ///
(clsn6 <- _cons@a, ) ///
(Slope@0 -> clsn1, ) /// // Determine temporal structure
(Slope@2 -> clsn3, ) ///
(Slope@3 -> clsn4, ) ///
(Slope@4 -> clsn5, ) ///
(Slope@5 -> clsn6, ) ///
, covstruct(_lexogenous, diagonal) ///
latent(Intercept Slope ) ///
cov(Slope@0 /// // Set slope variance to 0
e.clsn1@b ///
e.clsn3@b ///
e.clsn4@b ///
e.clsn5@b ///
e.clsn6@b) ///
means(Slope) /// // Estimate slope
nocapslatent
estimates store m3
estat gof, stats(chi2 rmsea indices residuals)
// Likelihood ratio test:
lrtest m3 m2 // Improvement in fit for Model 3
// Model 4: Random intercept, random slope (Table 2.6)
sem (Intercept@1 -> clsn1, ) /// // Constrain paths to be 1
(Intercept@1 -> clsn3, ) ///
(Intercept@1 -> clsn4, ) ///
(Intercept@1 -> clsn5, ) ///
(Intercept@1 -> clsn6, ) ///
(clsn1 <- _cons@a, ) /// // Constrain intercepts to be the same
(clsn3 <- _cons@a, ) ///
(clsn4 <- _cons@a, ) ///
(clsn5 <- _cons@a, ) ///
(clsn6 <- _cons@a, ) ///
(Slope@0 -> clsn1, ) /// // Determine temporal structure
(Slope@2 -> clsn3, ) ///
(Slope@3 -> clsn4, ) ///
(Slope@4 -> clsn5, ) ///
(Slope@5 -> clsn6, ) ///
, covstruct(_lexogenous, diagonal) ///
latent(Intercept Slope ) ///
cov(Intercept*Slope /// // Include intercept-slope covariance
e.clsn1@b ///
e.clsn3@b ///
e.clsn4@b ///
e.clsn5@b ///
e.clsn6@b) ///
means(Slope) /// // Estimate slope
nocapslatent
estimates store m4
estat gof, stats(chi2 rmsea indices residuals)
// Likelihood ratio test:
lrtest m4 m3 // Improvement in fit for Model 4
// Model 4 estimated as a multilevel model (Table 4.1):
// Read in data from clsn_cov.dat, this time using -corr2data-
#delimit ;
matrix input C = (
6.3944,
3.2716, 7.5282,
4.1435, 6.0804, 10.7290,
3.7058, 5.1597, 6.5672, 10.2920,
4.1286, 5.7608, 7.2365, 7.6463, 12.9085,
-0.0940, -0.0390, -0.1521, -0.1104, -0.1469, 0.2502
)
;
#delimit cr
corr2data clsn1 clsn3 clsn4 clsn5 clsn6 sex, n(851) ///
means(37.9542 37.2785 37.0463 36.5696 36.1363 0.49) ///
cov(C) cstorage(lower) clear
correlate, cov // Seems to work
gen id = _n // Create individual identifier
reshape long clsn, i(id) j(grade) // Convert to long format
mixed clsn grade || id: grade, var cov(uns)
Reference
Preacher, Kristopher J., Aaron L. Wichman, Robert C. MacCallum, and Nancy E. Briggs. 2008. Latent Growth Curve Modeling. Sage. doi: 10.4135/9781412984737
Labels:
#delimit,
corr2data,
estat gof,
estimates store,
Latent growth curve models,
lrtest,
matrix input,
mixed,
reshape,
sem,
ssd,
Textbooks
Aug 20, 2015
Random graphs (51): Bar graph
// Data from
// de Mouzon, Jacques, Veerle Goossens, Siladitya Bhattacharya, José Antonio
// Castilla, Anna Pia Ferraretti, V. Korsak, Markus Kupka, Karl-Gösta
// Nygren, Anders Nyboe Andersen, and the European IVF Monitoring
// Consortium for the European Society of Human Reproduction and Embryology.
// 2012. "Assisted Reproductive Technology in Europe, 2007. Results
// Generated from European Registers by ESHRE." Human Reproduction 27(4):
// 954-966. doi: 10.1093/humrep/des023
// Ferraretti, Anna Pia, Veerle Goossens, Jacques de Mouzon, Siladitya
// Bhattacharya, José Antonio Castilla, Vladislav Korsak, Markus Kupka, Karl
// G. Nygren, Anders Nyboe Andersen, and the European IVF Monitoring
// Consortium for the European Society of Human Reproduction and Embryology.
// 2012. "Assisted Reproductive Technology in Europe, 2008. Results
// Generated from European Registers by ESHRE." Human Reproduction 27(9):
// 2571-2584. doi: 10.1093/humrep/des255
// Ferraretti, Anna Pia, Veerle Goossens, Markus Kupka, Siladitya Bhattacharya,
// Jacques De Mouzon, José Antonio Castilla, Karin Erb, Vladislav Korsak,
// Anders Nyboe Andersen, and the European IVF Monitoring Consortium for the
// European Society of Human Reproduction and Embryology. 2013. "Assisted
// Reproductive Technology in Europe, 2009. Results Generated from European
// Registers by ESHRE." Human Reproduction 28(9):2318-2331. doi:
// 10.1093/humrep/det278
// Kupka, Markus S., Anna Pia Ferraretti, Jacques De Mouzon, Karin Erb, Thomas
// D'Hooghe, José Antonio Castilla, Carlos Calhaz-Jorge, Christian De
// Geyter, Veerle Goossens, and the European IVF Monitoring Consortium for
// the European Society of Human Reproduction and Embryology. 2014.
// "Assisted Reproductive Technology in Europe, 2010. Results Generated from
// European Registers by ESHRE." Human Reproduction 29(10):2099-2113. doi:
// 10.1093/humrep/deu175
// Always the two right-hands columns of Table 1
// "Cycles/million women 15–45 y" and "Cycles/million population"--Cycles of
// IVF, ICSI, and FER
// Fix country variable
kountry cntry, from(other) stuck marker
ren _ISO3N_ country
kountry country, from(iso3n) to(iso2c)
ren _ISO2C_ geo
drop country
ren cntry country
ren geo cntry
order cntry country
drop MARKER
// Fix string numbers
replace cycmio1545x2010 = subinstr(cycmio1545x2010, " ", "",.) // Remove blank (thousands separator)
replace cycmio1545x2009 = subinstr(cycmio1545x2009, " ", "",.)
replace cycmio1545x2008 = subinstr(cycmio1545x2008, " ", "",.)
replace cycmio1545x2007 = subinstr(cycmio1545x2007, " ", "",.)
destring cycmio1545x2010 cycmio1545x2009 cycmio1545x2008 cycmio1545x2007, replace
*list
// Turn data set around
reshape long cycmio1545x cycmio, i(cntry) j(year)
list, sepby(cntry)
// Cross-sectional bar graph at latest point in time
replace cycmio1545x = . if cntry == "DK" & year == 2010 // Get rid of dodgy DK 2010 value
by cntry: egen lastobs = max(cond(!missing(cycmio1545x), year, .)) // Identify year of latest observation
gen byte last = year == lastobs // ID variable of latest observation
// According to http://www.stata.com/support/faqs/data-management/first-and-last-occurrences/
graph hbar cycmio1545x if last == 1, over(country, sort(cycmio1545x) descending) nofill ///
ytitle("ART cycles per million women age 15{c 150}45", size(medsmall)) ylabel(0(2500)15000, ang(45)) ///
/*note("{it:Note:} Data refer to 2010 or latest available (no earlier than 2007)", span) */ ///
name(figure1, replace) ysize(6.5) xsize(5.5) scheme(lean1)
// Figure 1: ART cycles per million women age 15–45 per country, 2010
// Sources: Ferraretti et al. (2012; 2013) and Kupka et al. (2014).
// Notes: Values for Albania, Estonia, Ireland, Lithuania, Poland, Spain,
// Serbia, Switzerland, and Ukraine refer to 2008; for Croatia, Cyprus, France,
// and Denmark to 2009. ART cycles refer to IVF, ICSI, frozen embryo
// replacement (thawings), preimplantation genetic diagnosis and screening, egg
// donation (donation cycles), in vitro maturation, and frozen oocyte
// replacement (thawings).
Labels:
ESHRE,
graph hbar,
Identify last observation,
kountry,
Random graphs,
reshape
May 21, 2015
Random graphs (47): Complex scatterplot
import excel "C:\table+graph field match countries.xlsx", sheet("Sheet2") cellrange(A1:I27) clear firstrow
renvars, lower // Switch variable names to lower case
list // Check data
// Bring variable names into systematic shape
foreach x of varlist healthwelfare education socialsciencesbusinesslaw ///
artshumanities services sciencemathcomputing ///
agricultureveterinary engineeringmanufacturing {
ren `x' var_`x'
}
list // Check data
// Reshape data
reshape long var_, i(country) j(field_str) string
// Convert proportion into percentage
replace var_ = var_ * 100
// Get field_str variable into numerical format
encode field_str, gen(field1)
label define field1 1 "Agriculture and Veterinary" ///
2 "Arts and Humanities" ///
3 "Education" ///
4 "Engineering and Manufacturing" ///
5 "Health and Welfare" ///
6 "Science, Math, and Computing" ///
7 "Services" ///
8 "Social Sciences, Business, and Law", modify
// Sort by average segregation per field
list field1 var if country == "EU average"
recode field1 (5 = 1 "Health and Welfare") ///
(3 = 2 "Education") ///
(8 = 3 "Social Sciences, Business, and Law") ///
(2 = 4 "Arts and Humanities") ///
(7 = 5 "Services") ///
(6 = 6 "Science, Math, and Computing") ///
(1 = 7 "Agriculture and Veterinary") ///
(4 = 8 "Engineering and Manufacturing") ///
, gen(field)
// Plot
sort field var_
twoway (scatter field var_ if country != "EU average") ///
(scatter field var_ if country == "EU average", connect(l) lwidth(thick)) ///
, ylabel(1/8, val) xscale(alt) ///
xtitle("Percentage women by field of study" "in European countries") ///
ytitle("") note(" " "{it:Source:} European Labor Force Survey 2011, own calculations.", span) ///
legend(order(2) label(2 "Average" "across" "countries") ring(0) pos(1))
Jan 13, 2015
Random graphs (44): Small-multiples line plots
import excel "no of embryos transferred.xlsx", sheet("Tabelle1") firstrow clear
reshape long tot e1a e1p e2a e2p e3a e3p e4a e4p, i(country) j(year)
kountry country, from(other) stuck
ren _ISO3N_ ctry
kountry ctry, from(iso3n) to(iso2c)
drop ctry
ren _ISO2C_ cntry
replace cntry = "UK" if cntry == "GB"
order country cntry year tot e1a e2a e3a e4a e1p e2p e3p e4p
drop if cntry == "CY" | cntry == "EE" | cntry == "NL"
drop if cntry == "AL" | cntry == "BA" | cntry == "HR"
drop if cntry == "KZ" | cntry == "MD" | cntry == "RO"
drop if cntry == "UA" | cntry == "RS" | cntry == "TR"
drop if cntry == "ME"
twoway (line e1p year, by(cntry, note("")) cmissing(no) lwidth(thick)) ///
(line e2p e3p e4p year, by(cntry, note("{it:Sources:} ESHRE EIM reports, {it:Human Reproduction} 2001{c 150}2014", span)) cmissing(no)) ///
, xlabel(1997 2000 2005 2010, ang(h)) xtitle("") ytitle("% of all transfers") ///
legend(label(1 "1 embryo") label(2 "2 embryos") ///
label(3 "3 embryos") label(4 "4+ embryos") row(1))
Labels:
cmissing(no),
ESHRE,
import excel,
kountry,
Random graphs,
reshape,
twoway line
Dec 18, 2014
Random graphs (40): Line plot
clear
// Read in Table
input str4 cntry P1970 P1991 P1998 P2008 S1970 S1991 S1998 S2008 T1970 T1991 T1998 T2008
EU15 16.2 7.6 4.5 3.4 37.4 28.6 26.6 23.2 54 73 68.9 73.5
DE 8 3.4 2.5 2.1 48 35 30.4 25.3 44 62 67.1 72.5
FI 20 8 6.3 4.8 34 29 27.8 25.6 46 63 65.9 69.6
NL 6 4 3.6 3 36 24 19.9 16.7 58 72 76.5 80.3
SE 8 4 3.1 2.2 39 28 25.1 22.7 53 68 71.8 75.1
UK 3 2.1 1.9 1.5 42 25 . 17.7 55 73 75.3 80.7
CZ 13.5 10 5.6 3.5 49.4 45.9 41.4 38 37 44 53 58.6
end
// Crazy stuff to bring it into shape
preserve
drop S1970-T2008
reshape long P, i(cntry) j(year)
tempfile p
save `p', replace
restore
preserve
drop P1970-P2008 T1970-T2008
reshape long S, i(cntry) j(year)
tempfile s
save `s', replace
restore
drop P1970-S2008
reshape long T, i(cntry) j(year)
merge 1:1 cntry year using `p'
drop _merge
merge 1:1 cntry year using `s'
// Plot
twoway (scatter P year if cntry == "EU15", connect(L)) ///
(scatter S year if cntry == "EU15", connect(L)) ///
(scatter T year if cntry == "EU15", connect(L)) ///
, ytitle("Percentage of EU-15 workforce") ///
xtitle("") ///
note(" " "{it: Source:} Mau and Verwiebe 2010, pp. 153{c 150}4", span) ///
legend(label(1 "{bf:Primary sector:}" "Agriculture") ///
label(2 "{bf:Secondary sector:}" "Industry and" "construction") ///
label(3 "{bf:Tertiary sector:}" "Services")) ///
xlabel(1970 1991 1998 2008)
Labels:
input,
preserve,
Random graphs,
reshape,
tempfile,
twoway scatter
Nov 2, 2014
Random graphs (35): Line plots
clear
import excel "art births across countries.xlsx", sheet("Data") firstrow case(lower) clear
ren a geo // Rename country ID
// Rename variables: perc_YEAR format
foreach v of varlist b-o {
local x : variable label `v'
rename `v' perc`x'
}
// Reshape data set to long format
reshape long perc, i(geo) j(year)
// Cross-sectional bar graph at latest point in time
by geo: egen lastobs = max(cond(!missing(perc), year, .)) // Identify year of latest observation
gen byte last = year == lastobs // ID variable of latest observation
// According to http://www.stata.com/support/faqs/data-management/first-and-last-occurrences/
graph hbar perc if year >= 2007 & last == 1, over(geo, sort(perc) descending) nofill ///
ytitle("% ART live births among all births") ylabel(0(1)6) ///
note("{it:Note:} Data refer to 2010 or latest available (no earlier than 2007)", span) ///
/*title("Country differences in ART birth rates")*/ ///
name(crosssec, replace)
// Create classical twoway line plot
twoway (line perc year if geo == "DK", cmissing(no)) ///
(line perc year if geo == "SI", cmissing(no)) ///
(line perc year if geo == "BE", cmissing(no)) ///
(line perc year if geo == "SE", cmissing(no)) ///
(line perc year if geo == "DE", cmissing(no)) ///
(line perc year if geo == "FR", cmissing(no)) ///
(line perc year if geo == "IT", cmissing(no)) ///
(line perc year if geo == "UK", cmissing(no)) ///
(line perc year if geo == "US", cmissing(no)) ///
, xlabel(1997(1)2010, ang(h) labsize(small)) xtitle("") ///
ylabel(0(1)6) ytitle("% ART live births among all births") ///
legend(label(1 "DK") label(2 "SI") label(3 "BE") ///
label(4 "SE") label(5 "DE") label(6 "FR") ///
label(7 "IT") label(8 "UK") label(9 "US") ///
pos(2))
//cmissing(no) makes sure that gaps are not filled
// Use panel line plot functionality keep if inlist(geo, "DK", "SI", "BE", "SE", "DE", "FR", "IT", "UK", "US") encode geo, gen(country) xtset country year, yearly xtline perc, overlay xlabel(1997(1)2010, ang(h) labsize(small)) xtitle("") /// ylabel(0(1)6) ytitle("% ART live births among all births") ///
legend(pos(2)) // Note that the gaps are much harder to spot in this plot
// Small multiples version
twoway (line perc year, by(geo, note("") /*title("Development of ART birth rate, selected countries")*/) cmissing(no)), ///
xlabel(1998 (2) 2010, ang(h) labsize(small) alternate) xtitle("") ///
ylabel(0(1)6) yscale(r(0)) ytitle("% ART live births among all births") ///
name(longit, replace)
// yscale(r(0)) adds missing gridline according to
// http://www.stata.com/statalist/archive/2013-04/msg00904.html
Jul 2, 2014
Random graphs (24): Small-multiples and overlayed line plots
clear
// SG.GEN.PARL.ZS - Proportion of seats held by women in national parliaments (%)
wbopendata, indicator(SG.GEN.PARL.ZS) clear long
// Keep European countries
keep if regioncode == "ECS" // "Europe & Central Asia (all income levels)"
// Keep EU-28
keep if inlist(iso2code, "AT", "BE", "BG", "CY", "CZ", "DE", "DK", "EE", "ES") | ///
inlist(iso2code, "FI", "FR", "GB", "GR", "HR", "HU", "IE", "IT", "LI") | ///
inlist(iso2code, "LT", "LU", "LV", "NL", "PL", "PT", "RO", "SE", "SI") | ///
inlist(iso2code, "SK") // EU-28
// Identify and rename salient variables
ren sg_gen_parl_zs women
label var women "% of national parliament seats held by women"
ren iso2code cntry
drop countrycode region regioncode
// Reshape to drop empty observations
reshape wide women, i(cntry) j(year)
dropmiss women1960-women2014, force
reshape long
// Calculate EU-28 mean
preserve
collapse (mean) women, by(year)
gen cntry = "{bf:EU-28}"
tempfile eu28
save `eu28', replace
restore
append using `eu28'
// Plot small-multiples figure
twoway scatter women year ///
, by(cntry ///
, note(" ") ///
caption("{it:Source:} World Development Indicators, date of extraction: 2015-05-25", span size(small))) ///
cmissing(no) ///
connect(l) ///
msymbol(o) ///
xtitle("") xlabel(1990 1997 2000 2005 2010 2014, ang(45)) ///
ylabel(0 (10) 50) ///
name(small_multiples, replace)
// Plot overlayed figure
encode cntry, gen(country) // No strings for -xtline-
xtset country year, yearly // Declare panel
xtline women if cntry != "{bf:EU-28}" ///
, overlay ///
xlabel(1990 1997 2000 2005 2010 2014) xtitle("Year") ///
plot1opts(lpattern(dash)) plot2opts(lpattern(dash)) ///
plot3opts(lpattern(dash)) plot4opts(lpattern(dash)) ///
plot5opts(lpattern(dash)) plot6opts(lpattern(dash)) ///
plot7opts(lpattern(dash)) plot8opts(lpattern(dash)) ///
plot9opts(lpattern(dash)) plot10opts(lpattern(dash)) ///
plot11opts(lpattern(dash)) plot12opts(lpattern(dash)) ///
plot13opts(lpattern(dash)) plot14opts(lpattern(dash)) ///
plot15opts(lpattern(dash)) plot16opts(lpattern(dash)) ///
plot17opts(lpattern(dash)) plot18opts(lpattern(dash)) ///
plot19opts(lpattern(dash)) plot20opts(lpattern(dash)) ///
plot21opts(lpattern(dash)) plot22opts(lpattern(dash)) ///
plot23opts(lpattern(dash)) plot24opts(lpattern(dash)) ///
plot25opts(lpattern(dash)) plot26opts(lpattern(dash)) ///
plot27opts(lpattern(dash)) plot28opts(lpattern(dash)) ///
name(overlay, replace) ///
addplot(line women year if cntry == "{bf:EU-28}", ///
lwidth(thick) lpattern(solid) sort(cntry year)) ///
legend(order(29 "EU-28 (unweighted average)") pos(6) ring(0)) ///
caption("{it:Source:} World Development Indicators, date of extraction: 2015-05-25", span size(small)) ///
ytitle("% of national parliament seats" "held by women, EU-28 countries")
Labels:
append,
collapse,
dropmiss,
Random graphs,
reshape,
tempfile,
twoway scatter,
wbopendata,
xtline,
xtset
Subscribe to:
Posts (Atom)





























