/*
   This is quite a long Stata program file that provides a number of
     examples of visualizing data and statistics using the General
     Social Survey.

   Many commands, options, and sub-options can be used in different
     graphics.  These examples are not intended to be exhaustive,
     just illustrative.  You will need to read the online help file
     for more details.
*/

/* Read in a subset of the GSS data set */
#delimit ;
use year
    marital
    spkath
    spkrac
    spkhomo
    tvhours
    wordsum
    occ80
    abany
    sex
    race
    age
    educ
    degree
    marital
    prestige
    sexfreq
    partners
    nummen
    numwomen using http://terpconnect.umd.edu/~smilex3/GSS-Cumulative-72-12.dta, clear;
#delimit cr

/*
   Create two histograms, one bare-bones minimum,
     the other publication ready.
*/

/* Histogram #1--bare-bones */
hist marital, name(hist1, replace)

/* Histogram #2--publiction ready */
#delimit ;
hist marital, discrete percent gap(15)
  title("Distribution of Marital Status, 1972-2012", color(black))
  ytitle("Percent")
  xtitle("")
  note("Created for SOCY699C")
  caption("Source: General Social Survey")
  xlabel(1 "Married" 2 "Widowed" 3 "Separated" 4 "Divorced" 5 "NM", angle(0))
  ylabel(0(10)50, angle(0) nogrid)
  ymtick(0(5)50)
  color(black)
  fintensity(35)
  lwidth(vthin)
  lcolor(black)
  graphregion(color(white))
  name(hist2, replace);
#delimit cr


/*
   Bar garaph showing averages, which in this case is equal
     to percentages due to the recoding.
	
   Colors were determined using advice from colorbrewer
     (colorbrewer.org).  You can see the permanent link at:
      http://colorbrewer2.org/index.php?type=qualitative&scheme=Set1&n=3
*/
capture drop pctspkath pctspkrac pctspkhomo
recode spkath spkrac spkhomo (1=100) (2=0), pre(pct)

#delimit ;
graph bar pctspkath pctspkrac pctspkhomo if race<3 & year==2012, over(race)
  title("Percentage of Respondents Tolerant of Different Types" 
        "of People Speaking in Their Community", color(black))
  note("Source: General Social Survey, 2010")
  ylabel(, angle(0))
  bar(1, color(228  26  28))
  bar(2, color( 55 126 184))
  bar(3, color( 77 175  74))
  legend(label(1 "Anti-religionist")
         label(2 "Racist")
         label(3 "Homosexual")
         rows(1))
  graphregion(color(white))
  name(bar1, replace);
#delimit cr


/*
   Bar graph using a gray scale color scheme.  Note the
     colors -gs14-.  The -gs- series runs from -gs0- (black)
	 to -gs16- (white).  Also, the -glp- suboption (grid line
	 pattern) uses a dashed line for the y-grid.

   Also, note the use of the -gap(#)- suboption to the -over-
     option to control the distance between the bars.
*/
#delimit ;
graph bar tvhours if race<3, over(race, gap(20))
                             over(sex,  gap(80) relabel(1 "Men" 2 "Women"))
  title("Average Hours of TV Watched by Gender and Race", color(black))
  ytitle("Hours")
  note("Source: General Social Survey, 1975-2010")
  ylabel(, angle(0) glcolor(gs14) glp(dash))
  bar(1, color(gs14) lcolor(black))
  outergap(*.8)
  graphregion(color(white))
  name(bar2, replace);
#delimit cr


/*
   This example demonstrates how to manipulate the bar overlap,
     the legend, showing the values of each bar, and using the
     -asyvars- option which treats treats the first over() group
     as yvars.
*/
#delimit ;
graph bar wordsum if race<3, asy
    over(race)
    over(degree, relabel(1 "<HS" 2 "HS" 3 "JC" 4 "BA/BS" 5 "Grad"))
  title("Average Vocabulary Test Score")
  subtitle("by Race and Educational Degree")
  ytitle("Average # of correct words")
  note("Source: General Social Survey, 1972-2010")
  ylabel(, angle(0))
  legend(rows(1) ring(0) position(11))
  bargap(-30)
  blabel(bar, position(inside) format(%9.1f) color(white))
  name(bar3, replace);
#delimit cr


/*
   If you have a lot of categories (e.g. year), consider
     using a dot chart.
*/
#delimit ;
graph dot tvhours if race<3, asyvars nofill
    over(race)
    over(year, label(labsize(small)))
  title("Average Hours of Television Watched in a Day")
  ytitle("Hours")
  note("Source: General Social Survey")
  exclude0
  ylabel(2.5(.5)4.5, format(%3.1f))
  ymtick(2.5(.1)4.5)
  legend(col(1) ring(0) position(2))
  name(dot1, replace);
#delimit cr


/*
   OCC80
   managerial and professional specialty occupations
   technical, sales, and administrative support occupations
   service occupations
   precision production, craft, and repair occupations
   operators, fabricators, and laborers
*/
capture drop occupation
#delimit ;
recode occ80 (  3/199=1 "Managers")
             (203/389=2 "Technical")
             (403/499=3 "Service")
             (503/699=4 "Precision")
             (703/899=5 "Laborers")
             (else=.), gen(occupation);
#delimit cr


/* Another good way to compare the distribution of a variable
     over discrete categories is to use box-and-whisker or
     box plot.  There is a good description at
     http://en.wikipedia.org/wiki/Box_plot.
*/
graph box wordsum, over(occupation) name(box1, replace)

#delimit ;
graph box wordsum if race<3, asy nooutsides
    over(race)
    over(occupation, axis(noline))
  title("Vocabulary Test Score Distribution", color(black))
  subtitle("Comparing Race within Occupation")
  ytitle("# of correct words")
  note("Source: General Social Survey, 1988-2010")
  ylabel(, angle(0) glpattern(dash))
  yscale(off)
  legend(region(color(white)))
  graphregion(color(white))
  name(box2, replace);
#delimit cr

/*
   This example uses the -collapse- command to create create
     two series showing the proportion of white and black
     respondents who favor legal abortion.

   I have suppressed the legend but have used the -text- option
     to differentiate the two series.  You can use this option to
     place text anywhere on your graphic.
*/
capture drop abanydum
gen abanydum=2-abany

capture drop racedum
recode race (1=1 "white") (2/3=0 "non-white"), gen(racedum)

preserve
  collapse abanydum, by(year racedum)
  drop if abanydum==.
  list
  #delimit ;
  graph twoway (line abanydum year if racedum==0, lw(thick))
               (line abanydum year if racedum==1, lw(thick)),
  title("Proportion of People Favoring Legal Abortions by Race")
  xtitle("") ytitle("")
  note("Source: General Social Survey, 1977-2010")
  legend(off)
  xlabel(1980(5)2010) ylabel(, format(%5.2f))
  text(0.35 2006 "Black") text(0.45 2006 "White");
restore


/*
   This example shows how you can combine graphs and use
     different axes for scaling, one on the left and one
     on the right (also, top and bottom if you want)

   Note the use of the -separate- command which often comes
     in handy.  I need to use separate here so I can calculate
     the black/white difference, unlike the previous exmaple.

   There is a lot of options being used so be sure to read the
     onlinehelp file for all o fthe details.
*/
preserve
  separate abanydum, by(racedum)
  collapse abanydum?, by(year)
  drop if missing(abanydum0, abanydum1)

  gen diff=abanydum1-abanydum0
  list

  #delimit ;
  graph twoway (line abanydum0 year, lw(thick) yaxis(1 2) xaxis(1 2))
               (line abanydum1 year, lw(thick) yaxis(1 2) xaxis(1 2))
               (line diff      year, lw(thick) lc(midblue) yaxis(1 2))
               (lfit diff      year, yaxis(1 2) lc(black) lp(dash)),
    title("Change in Proportion Favoring Legal"
       "Abortion for Any Reason by Race")
        xtitle("") xtitle("", axis(2))
        ytitle("Proportion")
        note("Source: General Social Survey")
        xlabel(1980(10)2010, axis(1))
    xlabel(1980(10)2010, axis(2) nolabels noticks)
    ylabel(0(.05).1 , axis(2) gmin angle(horizontal) format(%4.2f) grid)
    ylabel(.3(.05).5, axis(1) gmax angle(horizontal) format(%4.2f))
    legend(label(1 "Black")
           label(2 "White")
           label(3 "Difference")
           label(4 "OLS") 
			  col(1) pos(3) size(small) region(lstyle(none)))
    graphregion(color(white))
    name(twoscalegraph, replace);
  #delimit cr
restore


/*
   Sometimes it is conceptually easier to create a new grouping
     variable instead of using multiple -over- options.  The -egen-
     command has many forms and is well worth exmaining in detail.

   I use the -group- form which creates a variable taking on
     values 1, 2, ... for the groups formed by varlist.  In this
     case this produces the intersection of race and sex, 1=non-white male,
     2=non-white female, 3=white male, and 4=white female.

   This mesure can now easily be used for categories in graphs.
*/
capture drop rs
egen rs=group(racedum sex), label
tab rs

preserve
  collapse abanydum, by(rs)
  list
  graph bar abanydum, over(rs)
restore


replace nummen  =. if nummen   > 989
replace numwomen=. if numwomen > 989
gen partners1=nummen+numwomen
label variable partners1 "Number of sexual partners since age 18"

preserve
  collapse partners1, by(year rs)
  drop if partners1==.
  list
  #delimit ;
  graph twoway (line partners1 year if rs==1, lw(thick))
               (line partners1 year if rs==2, lw(thick))
               (line partners1 year if rs==3, lw(thick))
               (line partners1 year if rs==4, lw(thick)),
    title("Number of sexual partners since age 18")
        subtitle("Comparing race and gender")
        xtitle("") ytitle("# of partners")
        ylabel(, angle(0)) ymtick(0(25)30)
        legend(region(lcolor(white)) row(1)
           label(1 "nw men")
           label(2 "nw women")
           label(3 "w men")
           label(4 "w women"))
    graphregion(color(white))
    name(partners1, replace);
restore


/*
   Some graphics related to regression analyses.  Here are linear and quadratic
     fitted lines with 95% confidence intervals for the predicted values.

   The first example creates two graphs but suppresses their output.  Then,
     they are combined into a single graph to allow better comparison.
*/
capture drop sexfreq1
recode sexfreq (0=0) (1=2) (2=12) (3=36) (4=52) (5=156) (6=208), gen(sexfreq1)
regress sexfreq1 age

#delimit ;
graph twoway (lfitci sexfreq1 age if year==2010)
             (lfit   sexfreq1 age if year==2010),
  title("Linear Fit")
  ylabel(0(25)100)
  yscale(range(0 100))
  legend(off)
  aspectratio(1)
  nodraw
  name(regress1, replace);

graph twoway (qfitci sexfreq1 age if year==2010)
             (qfit   sexfreq1 age if year==2010),
  title("Quadratic Fit")
  ylabel(0(25)100)
  yscale(range(0 100))
  legend(off)
  aspectratio(1)
  nodraw
  name(regress2, replace);

graph combine regress1 regress2,  nocopies;
#delimit cr


/*
   Estimate a multiple regresison model and show a number
     of plots of predicted values.
*/
regress sexfreq1 c.age##c.age i.marital if year==2010
margins, over(marital)
marginsplot, recast(scatter) name(margins1, replace)

regress sexfreq1 c.age##c.age i.marital
margins, over(year)
marginsplot, recastci(rarea) recast(line) name(margins2, replace)

margins, over(year marital)
marginsplot, recastci(rarea) recast(line) name(margins3, replace)


/*
   Margins can be used to estimate the expected probabilities
     in a logistic regression model.
*/
logistic abanydum racedum educ age i.marital
margins, over(year racedum)
marginsplot, recast(line) recastci(rarea) xlabel(1980(10)2010) name(margins4, replace)