*! xtgraph.ado, Version 0.31 *! by PT Seed (paul.seed@kcl.ac.uk) *! graphs of xt (cross-sectional time-series) data *! Presents averages with error bars of a single variable, by time *! data can be grouped. * 14/4/2001 All major functions now seem to be working, ready for SUG * Log & power transformations, model fitting, outputting data & graphs, * The only problem seems to be that the undeclared "show" does not * work with axis labels on the graph. * -av(am)- now allowed as a synonym for -av(mean)- prog define xtgraph set trace off version 6.0 *di "***** parse the command *************** error 0 syntax varlist(min=1 max=1 numeric) [if] [in] , /* */ [AVerage(string) bar(string) /* */ POWer(real 1) log(string) model /* */ i(string) t(string) GRoup(varlist max=1) Half /* */ OFFset(real 0) show MINobs(real 1) level(real $S_level) /* */ list savedat(string) by(string) listwise missing /* */ L1title(string) Symbol(string) Connect(string) Weight(string) nograph saving(string) * ] *di "***** ***** ***** ***** ***** ***** *di "***** Collect the I and T vars ***** xt_iis `i' local ivar "`s(ivar)'" xt_tis `t' local tvar "`s(timevar)'" *di "***** ***** ***** ***** ***** ***** *di "***** Handle the conditions ***** preserve local yvar "`varlist'" cap keep `if' `in' qui keep if `yvar' ~= . cap keep if `group' ~= . tempvar n if `minobs' > 1 { qui egen `n' = count(`yvar'), by(`group' `by' `tvar') cap keep if `n' >= `minobs' if _rc & "`1'" ~= "" { di in red "Invalid option minobs(" in ye"`minobs'" in red")" exit 198 } drop `n' } if "`listwise'" ~= "" { qui egen `n' = count(`yvar'), by(`ivar') summ `n' , meanonly keep if `n' == r(max) drop `n' } if "`missing'" ~= "" { drop if `group' == . if "`by'" ~= "" {drop if `by' == .} } qui egen `n' = count(`yvar'), by(`group' `by' `tvar') if "`saving'" ~= "" { local graph = "" } ***** ***** ***** ***** ***** ***** *di "***** Check the syntax ***** if "`log'" ~= "" { cap confirm number `log' if _rc == 7 { di in red "Option log() must contain a number." exit 198 } } if "`log'" ~= "" & "`power'" ~= "1" { di in red "Options log() and power() canot be used together exit 198 } if "`average'" == "" | "`average'" == "am" {local average = "mean" } if ("`average'" == "gm" | "`average'" == "hm" ) & "`power'" ~= "1" { di in red "Option power() should not be used with average(`average')" exit 198 } if "`average'" == "gm" { local average "mean" local power = 0 local log = 0 } if "`average'" == "hm" { local average "mean" local power = -1 } if "`bar'" == "" {local bar "ci" } if "`average'" == "median" { if "`bar'" ~= "ci" & "`bar'" ~= "no" & "`bar'"~= "iqr" & "`bar'" ~= "rr" { di in red "Option bar(`bar') invalid with average(`average')." exit 198 } } if "`average'" == "mean" { if "`bar'" ~= "ci" & "`bar'" ~= "no" & "`bar'"~= "se" /* */ & "`bar'"~= "sd" & "`bar'" ~= "rr" { di in red "Option bar(`bar') invalid with average(`average')." exit 198 } } if "`group'" == "" { tempvar group qui gen `group' = 1 } *di "***** Symbol ***** tempvar gpno sort `group' qui by `group' : gen `gpno' = _n == 1 qui replace `gpno' = sum(`gpno') sort `gpno' local ngrp = `gpno'[_N] while length("`symbol'") < `ngrp' { local symbol "`symbol'OSTodp" } local symbol = substr("`symbol'",1,`ngrp') local symbol "`symbol'ii" *di "***** Connect ***** while length("`connect'") < `ngrp' { local connect "`connect'l" } local connect "`connect'II" *di "***** Half ***** if "`half'" ~= "" & `ngrp' ~= 2 { di in red "Half-bars are not appropriate with more than two groups." di in red "Perhaps you should use the offset()option." exit 198 } *di "***** Collect the Y and T labels ***** local y_vrlab : var label `varlist' if "`y_vrlab'" == "" { local y_vrlab = "`varlist'" } local y_vllab : val label `varlist' local t_vrlab : var label `tvar' if "`t_vrlab'" == "" { local t_vrlab = "`tvar'" } local t_vllab : val label `tvar' ***** ***** ***** ***** ***** ***** * di "***** Carry out the transformations ***** if "`power'" ~= "1" & "`power'" ~= "0" & "`model'" == "" { qui replace `yvar' = `yvar'^`power' } else if "`log'" ~= "" & "`model'" == "" { cap assert `yvar' - `log' > 0 if _rc { local ln_sign = "-" cap assert -`yvar' - `log' > 0 if _rc { di in red "Non-positive values of `yvar' - `log' encountered." exit _rc } } qui replace `yvar' = log(`ln_sign'1*`yvar'-`log') } ****************************************** *di "***** Calculate the averages and bars *********** tempvar av if "`bar'" ~= "no" { tempvar lb ub } *di " ***** ***** ***** means ***** ***** ***** if "`average'" == "mean" { if "`model'" == "" { qui egen `av' = `average'(`yvar'), by(`group' `by' `tvar') } else { cap predict `av' if _rc { di in red "The last model fitted cannot be detected." exit _rc } } *di " ***** ***** sd & rr ***** ***** if "`bar'" == "sd" | "`bar'" == "rr" { tempvar sd if "`model'" == "" { qui egen `sd' = sd(`yvar'), by(`group' `by' `tvar') } else { cap predict `sd' , stdf if _rc { di in red "The last model fitted does not permit estimates of stdp, di in red "needed for standard deviation bars" exit _rc } } if "`bar'" == "sd" { qui gen `lb' = `av' - `sd' qui gen `ub' = `av' + `sd' } if "`bar'" == "rr" { qui gen `lb' = `av' - invnorm(.5+`level'/200)*`sd' qui gen `ub' = `av' + invnorm(.5+`level'/200)*`sd' } } *di " ***** ***** se & ci ***** ***** if "`bar'" == "ci" | "`bar'" == "se" { tempvar se if "`model'" == "" { qui egen `se' = sd(`yvar'), by(`group' `by' `tvar') qui replace `se' = `se'/`n'^.5 } else { qui predict `se', stdp if _rc { di in red "The last model fitted does not permit estimates of stdf, di in red "(the standard deviaiton of the forecast) needed for standard error bars" exit _rc } } if "`bar'" == "ci" { qui gen `lb' = `av' - (invt(`n'-1,`level'/100)*`se') qui gen `ub' = `av' + (invt(`n'-1,`level'/100)*`se') } if "`bar'" == "se" { qui gen `lb' = `av' - `se' qui gen `ub' = `av' + `se' } } } ***** ***** ***** ***** ***** ***** ***** *di "***** ***** ***** medians ***** ***** ***** if "`average'" == "median" & "`bar'" == "ci" { qui gen `av' = . qui gen `lb' = . qui gen `ub' = . sort `group' `by' `tvar' tempvar gp_t qui by `group' `by' `tvar': gen `gp_t' = _n == 1 qui replace `gp_t' = sum(`gp_t') local n_gp_t = `gp_t'[_N] local i = 1 while `i' <= `n_gp_t' { qui centile `yvar' if `gp_t' == `i', level(`level') qui replace `av' = r(c_1) if `gp_t' == `i' qui replace `lb' = r(lb_1) if `gp_t' == `i' qui replace `ub' = r(ub_1) if `gp_t' == `i' local i = `i' + 1 } } if "`average'" == "median" & "`bar'" == "iqr" { qui egen `av' = pctile(`yvar') , by(`group' `by' `tvar') qui egen `lb' = pctile(`yvar') , by(`group' `by' `tvar') p(25) qui egen `ub' = pctile(`yvar') , by(`group' `by' `tvar') p(75) } if "`average'" == "median" & "`bar'" == "rr" { local l_cent = 50-`level'/2 local u_cent = 50+`level'/2 qui egen `av' = pctile(`yvar') , by(`group' `by' `tvar') qui egen `lb' = pctile(`yvar') , by(`group' `by' `tvar') p(`l_cent') qui egen `ub' = pctile(`yvar') , by(`group' `by' `tvar') p(`u_cent') } ***** ***** ***** ***** ***** ***** * di "***** Carry out the back transformations ***** if "`power'" ~= "1" & "`power'" ~= "0" { qui replace `av' = `av'^(1/`power') if "`bar'" ~= "no" & `power' >=0 { qui replace `lb' = `lb'^(1/`power') qui replace `ub' = `ub'^(1/`power') } if "`bar'" ~= "no" & `power' <0 { tempvar lb_temp qui gen `lb_temp' = `ub'^(1/`power') qui replace `ub' = `lb'^(1/`power') qui replace `lb' = `lb_temp' } } if "`log'" ~= "" { qui replace `av' = `ln_sign'1*(exp(`av')+`log') if "`bar'" ~= "no" { qui replace `lb' = `ln_sign'1*(exp(`lb')+`log') qui replace `ub' = `ln_sign'1*(exp(`ub')+`log') } } if "`power'" == "0" & "`log'" == "" { qui replace `av' = exp(`av') qui replace `lb' = exp(`lb') qui replace `ub' = exp(`ub') } ***** ***** ***** ***** ***** *di "***** Offset the tvar ***** tempvar t2 if `offset' ~= 0 { qui tab `gpno', qui gen `t2' = `tvar' + (`gpno'-(r(r)+1)/2)*`offset' } else { qui gen `t2' = `tvar' } label var `t2' "`t_vrlab'" label val `t2' `t_vllab' ***** ***** ***** ***** ***** *di "***** Drop the extra variables & values ***** sort `group' `by' `tvar' tempvar touse qui by `group' `by' `tvar' : gen `touse' = _n == 1 qui keep if `touse' keep `ivar' `tvar' `n' `group' `gpno' `by' `yvar' `av' `avlist' `lb' `ub' `t2' if `group'[_N] == `group'[1] { drop `group' local group } ***** ***** ***** ***** ***** *di "***** Rename the variables ***** cap rename `lb' lb if _rc == 0 { local lb lb } cap rename `ub' ub if _rc == 0 { local ub ub } cap rename `t2' t2 if _rc == 0 { local t2 t2 } cap rename `av' `average' if _rc == 0 { local av `average'} cap noi rename `n' n if _rc == 0 { local n n} ***** ***** ***** ***** ***** *di "***** Label the groups ***** cap assert "`group'" ~= "" if _rc { local avlist "`av'" label var `avlist' "`y_vrlab'" local c "l" } else { tempvar grp_num cap confirm numeric var `group' if _rc == 0 { local grp_num "`group'" } else { encode `group' , gen(`grp_num') } * grp_num is the numeric form of group. local i = 1 while `i' <= `ngrp' { qui summ `grp_num' if `gpno' == `i', meanonly local gr_val = r(mean) local label : label (`grp_num') `gr_val' * First try to make name for average from label local name = lower(substr("`label'",1,8)) while index("`name'", " ") { local l_name = substr("`name'", 1, index("`name'", " ")-1) local r_name = substr("`name'", index("`name'", " ")+1, .) local name "`l_name'_`r_name'" } cap confirm new var `name' while _rc { * Then try adding underscores local j = `j' + 1 local name = "_" + substr("`name'",1,7) cap confirm new var `name' if `j' == 8 { * Then try numbered groups di in ye "Problem giving meaningful name to group `i'." di in ye "Average is numbered only." local name "Group_`gr_val'" } if `j' == 10 { * Then give up di in red "Problem giving meaningful names to all groups." di "No solution." exit 198 } *End of [not] new var `name' } *di "***** Rename the group average for group `i' ***** qui gen `name' = `av' if `gpno' == `i' label var `name' "`label'" local avlist "`avlist' `name'" local c "`c'l" local i = `i' + 1 * End of consideration of name for group `i' } * End of consideration of group names for multiple groups } ***** ***** ***** ***** ***** ***** *di "***** Label the Y axis ***** if "`l1title'" == "" { local l1title "`y_vrlab'" } local l1title "l1(`l1title')" *di "***** List the data ***** if "`list'" ~= "" { if "`average'" == "mean" { if `power' == 0 | "`log'" == "0" { local av_type "Geometric mean" } else if `power' == -1 { local av_type "Harmonic mean" } else if "`power'" ~= "1" { local av_type = "Transformed mean [x^`power']" } else if "`log'" ~= "" { local av_type = "Transformend mean [log(x -`log')]" } else { local av_type "Arithmetic mean" } } else local av_type = "`average'" di _n "`av_type' of `varlist' " _c if "`bar'" == "ci" | "`bar'" == "rr" { local l = "`level'% " } if "`bar'" ~= "no" { local up_bar = upper("`bar'") di "with bars based on `l'`up_bar'" } if "`group'" ~= "" & "`by'" == "" { di "by `group' and `tvar'." } else if "`group'" ~= "" & "`by'" ~= "" { di "by `group', `by' and `tvar'." } else if "`group'" == "" & "`by'" ~= "" { di "by `by' and `tvar'." } else di "by `tvar'." if "`group'" ~= "" | "`by'" ~= "" { sort `group' `by' by `group' `by' : list `tvar' `n' `av' `lb' `ub' , noobs } else { list `tvar' `n' `av' `lb' `ub' , noobs } } ***** ***** ***** ***** ***** *di "***** Weight the points ***** if "`weight'" ~= "" { if "`weight'" === "n" { local wtvar "`n'" } else if "`weight'" == "bar" { tempvar wtvar qui gen wtvar = `ub' - `lb' } else local wtvar = `weight' local weight "[fw=`wtvar']" } *di "***** Save the data ***** if "`savedat'" ~= "" { tempfile temp qui save `temp' keep `group' `by' `tvar' `n' `av' `lb' `ub' `wtvar' save `savedat' qui use `temp', replace } ***** ***** ***** ***** ***** *di "***** Handle the half-bars ***** if "`half'" ~= "" { tempvar lbs ubs qui gen str1 `lbs' = "" qui gen str1 `ubs' = "" sort `by' `tvar' `av' `group' qui by `by' `tvar' : replace `lb' = `av' if `group' == `group'[_N] qui by `by' `tvar' : replace `ub' = `av' if `group' == `group'[1] } ***** ***** ***** ***** ***** *di "***** Sort the data ***** sort `by' `group' `t2' ***** ***** ***** ***** ***** *di "***** Draw the graph ***** if "`by'" ~= "" { local by "by(`by')" } if "`saving'" ~= "" { local saving "saving(`saving')" } if "`show'" ~= "" { cap noi di "graph `avlist' `lb' `ub' `t2' `weight', connect(`connect') symbol(`symbol') `xlab' `ylab' `by' `options' `saving' " } if "`graph'" == "" {graph `avlist' `lb' `ub' `t2' `weight', connect(`connect') symbol(`symbol') `xlab' `ylab' `l1title' `by' `saving' `options'} ************************************* end xtgraph ************************************* ************************************* ************************************* *! version 2.2.0 01jul1999 program define centile, rclass version 6 syntax [varlist] [if] [in] [, CCi /* */ Centile(numlist >=0 <=100) /* */ Level(real $S_level) Meansd Normal ] tempvar touse notuse mark `touse' `if' `in' qui gen byte `notuse' = . /* will be reset for each variable */ if "`centile'"=="" { local centile 50 } /* Parse `centile' and count the requested centiles, placing them into c1, c2, .... */ local nc 0 tokenize "`centile'" while "`1'" != "" { local nc = `nc' + 1 local c`nc' `1' local cents "`cents' `1'" /* to return in r() */ mac shift } local tl1 " Obs " local ttl "Percentile" if "`meansd'"=="" { if "`normal'"=="" { if "`cci'"~="" { di in gr _n _col(52) "-- Binomial Exact --" } else { di in gr _n _col(52) "-- Binom. Interp. --" } } else { di in gr _n _col(32) /* */ "-- Normal, based on observed centiles --" } } else { di in gr _n _col(32) "-- Normal, based on mean and std. dev.--" } #delimit ; di in gr "Variable | `tl1' `ttl' Centile [`level'% Conf. Interval]" _n "---------+" _dup(61) "-" ; #delimit cr local anymark 0 local alpha2 = (100-`level')/200 local zalpha2 = -invnorm(`alpha2') /* Run through varlist */ tokenize `varlist' local vl while "`1'" != "" { capt conf str var `1' if _rc { local vl "`vl' `1'" } mac shift } local nobs 0 /* in case loop aborted */ tokenize `vl' while "`1'" ~= "" { local yvar "`1'" /* notuse: 0 = use, 1 = not use -- sort will put use first */ qui replace `notuse' = ~`touse' qui replace `notuse' = 1 if `yvar'==. sort `notuse' `yvar' qui sum `yvar' if ~`notuse' local nobs = r(N) local mean = r(mean) local sd = sqrt(r(Var)) /* Formatting */ local skip = 8 - length("`yvar'") local fmt : format `yvar' if substr("`fmt'",-1,1)=="f" { local ofmt="%9."+substr("`fmt'",-2,2) } else if substr("`fmt'",-2,2)=="fc" { local ofmt = "%9." + substr("`fmt'",-3,3) } else local ofmt "%9.0g" /* Calc required centile(s) */ local j 1 local s 7 while `j' <= `nc' { local mark "" local cj "c`j'" local quant = ``cj''/100 if "`meansd'" ~= "" & (`nobs' > 0) { /* Normal distribution (parametric estimates) */ local z = invnorm(`quant') local centil = `mean'+`z'*`sd' local se = `sd'*sqrt(1/`nobs'+(`z')^2/(2*(`nobs'-1))) local cLOWER = `centil'-`zalpha2'*`se' local cUPPER = `centil'+`zalpha2'*`se' } else if `nobs' > 0 { /* If `normal' is not set, calc centile and exact nonparametric confidence limits (for example, see Gardner and Altman 1989, pp 72-74), interpolating when not at ends of distribution. An iterative process is required for each limit. As starting values, find ranks for lower and upper limits using a normal approximation. If `normal' is set, use normal distribution formula for variance, (10.29) in Kendall & Stuart (1969) p 237. `alpha2' is for example .025 for a 95% CI. */ local frac1 = (`nobs'+1)*`quant' * di "Central fraction and rank = " `quant',`frac1' local i1 = int(`frac1') local frac1 = `frac1'-`i1' if `i1' >= `nobs' { local centil = `yvar'[`nobs'] } else if `i1' < 1 { local centil = `yvar'[1] } else { local centil = `yvar'[`i1']+ /* */ `frac1'*(`yvar'[`i1'+1]-`yvar'[`i1']) } if "`normal'" == "" { local nq = `nobs'*`quant' local z = sqrt(`nq'*(1-`quant'))*`zalpha2' local rzLOW = int(.5+`nq'-`z') local rzHIGH = 1+int(.5+`nq'+`z') * di "lower and upper approx ranks = " `rzLOW',`rzHIGH' local r1 `rzHIGH' if `r1' > `nobs'+1 { local r1 = `nobs'+1 } local r0 = `r1'-1 local p0 = Binomial(`nobs',`r0',`quant') local p1 = Binomial(`nobs',`r1',`quant') local done 0 while ~`done' { if `p0' > `alpha2' { if `p1' <= `alpha2' { local done 1 } else { local r0 = `r1' local p0 = `p1' local r1 = `r1'+1 local p1 = Binomial(`nobs',`r1',`quant') } } else if `p0' == `alpha2' { local r1 = `r0' local p1 = `p0' local done 1 } else { local r1 = `r0' local p1 = `p0' local r0 = `r0'-1 local p0 = Binomial(`nobs',`r0',`quant') } } * di "Upper r0, p0, r1, p1, interp =" `r0',`p0',`r1',`p1',(`p0'-`alpha2')/(`p0'-`p1') /* Upper limit. Interpolate between r0 and r1, r1 being conservative. Note that p0>=p1 (upper tail areas). */ if `r0' >= `nobs' { local cUPPER = `yvar'[`nobs'] local mark "*" local anymark 1 } else if `r0' < 1 { local cUPPER = `yvar'[1] local mark "*" local anymark 1 } else { if "`cci'" == "" { local cUPPER = `yvar'[`r0'] /* */ +((`p0'-`alpha2')/(`p0'-`p1')) /* */ *(`yvar'[`r1']-`yvar'[`r0']) } else { local cUPPER = `yvar'[`r1'] } } local r1 `rzLOW' if `r1' < 0 { local r1 0 } local r0 = `r1'-1 local p0 = 1-Binomial(`nobs',`r0'+1,`quant') local p1 = 1-Binomial(`nobs',`r1'+1,`quant') local done 0 while ~`done' { if `p1' > `alpha2' { if `p0' <= `alpha2' { local done 1 } else { local r1 = `r0' local p1 = `p0' local r0 = `r0'-1 local p0 = 1-Binomial(`nobs',`r0'+1,`quant') } } else if `p0' == `alpha2' { local r0 = `r1' local p0 = `p1' local done 1 } else { local r0 = `r1' local p0 = `p1' local r1 = `r1'+1 local p1 = 1-Binomial(`nobs',`r1'+1,`quant') } } /* Interpolate between r1 and r1+1, r1 being conservative. Note that p0<=p1 (both are lower tail areas). */ if `r1' >= `nobs' { local cLOWER = `yvar'[`nobs'] local mark "*" local anymark 1 } else if `r1' < 1 { local cLOWER = `yvar'[1] local mark "*" local anymark 1 } else { if "`cci'" == "" { local cLOWER = `yvar'[`r1'] /* */ +((`alpha2'-`p0')/(`p1'-`p0')) /* */ *(`yvar'[`r1'+1]-`yvar'[`r1']) } else { local cLOWER = `yvar'[`r1'] } } * di "Lower r0, p0, r1, p1, interp =" `r0',`p0',`r1',`p1',(`alpha2'-`p0')/(`p1'-`p0') * di "Var = `yvar', centile `cj'=" `centil' ", CI = " `cLOWER',`cUPPER' } else { /* Normal distribution, observed centiles */ local dens = exp(-0.5*((`centil'-`mean')/`sd')^2)/(`sd'*sqrt(2*_pi)) local se = sqrt(`quant'*(1-`quant')/`nobs')/`dens' local cLOWER = `centil'-`zalpha2'*`se' local cUPPER = `centil'+`zalpha2'*`se' } } if (`j' == 1) & (`nobs' > 0) { di in gr _skip(`skip') "`yvar' |" _col(10) /* */ in yel %8.0f `nobs' /* */ _col(23) %7.0g ``cj'' /* */ _col(35) `ofmt' `centil' /* */ _col(51) `ofmt' `cLOWER' /* */ _col(63) `ofmt' `cUPPER' in gr "`mark'" } else if `nobs' > 0 { di in gr " |" /* */ in yel _col(23) %7.0g ``cj'' /* */ _col(35) `ofmt' `centil' /* */ _col(51) `ofmt' `cLOWER' /* */ _col(63) `ofmt' `cUPPER' in gr "`mark'" } else { /* `nobs' == 0 */ if (`j' == 1) { di in gr _skip(`skip') "`yvar' |" /* */ _col(10) in yel %8.0f `nobs' } local centil . local cLOWER . local cUPPER . } /* Store centile in S_# macros, starting at 7. Store centiles also in r(c_#) where # starts at 1. Store lower and upper bounds on centiles in r(lb_#) and r(ub_#) */ local tmp = `s' - 6 ret scalar c_`tmp' = `centil' /* save in r() */ global S_`s' `centil' /* also save in S_# */ /* save confidence limits also in r(), but not S_# */ ret scalar lb_`tmp' = `cLOWER' ret scalar ub_`tmp' = `cUPPER' local j = `j'+1 local s = `s'+1 } mac shift } if "`anymark'" == "1" { di in gr _n /* */ "`mark' Lower (upper) confidence limit held at minimum (maximum) of sample" } /* Store quantities at final point; S_4 is duplicate of S_[`nc'+6]. */ ret scalar N = `nobs' ret scalar n_cent = `nc' ret local centiles `cents' /* double save in S_# */ global S_1 `nobs' global S_2 `nc' global S_3 ``cj'' global S_4 `centil' global S_5 `cLOWER' global S_6 `cUPPER' end exit