-------------
      name:  <unnamed>
       log:  /Users/eric/Work/MB/ATCR/200/Labs/Lab 1 exercises 2018.smcl
  log type:  smcl
 opened on:   5 Jul 2018, 05:33:02

.
. /* Read in the lab1 data posted on the class web site.  These data come from a
> study of voluntary counseling and testing (VCT) for HIV at Mulago Hospital
> in Uganda. */
.
. use lab1.dta, clear

.
. /* We will do a little data management to make the data easier to use. First
> assign some variable labels */
.
. label var age "Age"

. label var cd4 "CD4 count"

. label var cd4_cat "CD4 count"

. label var hiv "HIV status"

. label var hungry "How often household members go hungry"

. label var lastalc "Last time used alcohol"

. label var sex "Sex"

. label var n_adults "# adults in household"

. label var n_child "# people <18yo in household"

. label var n_support "# people supported financially"

. label var religion "Religion"

.
. /* Now add some missing value labels. A easy-to-remember convention is to
> use the variable name as the name of the value label. (Stata doesn't have
> any problem with this). */
.
. label define sex 1 "Male" 2 "Female"

. label values sex sex

. label define hiv 0 "HIV-" 1 "HIV+"

. label values hiv hiv

.
. * Stata lets us look at a spreadsheet of the data
.
. browse

.
. /* Exploratory data analysis -- descriptive statistics.
>
> First see what variables are in the dataset.  Check how many observations
> and variables the dataset includes. */
.
. describe

Contains data from lab1.dta
  obs:         3,389
 vars:            12                          4 Jul 2018 08:10
 size:        44,057
-------------
              storage   display    value
variable name   type    format     label      variable label
-------------


age               byte      %12.0g                      Age


auditc            byte      %12.0g                      AUDITC, past 3 mos


cd4               int       %11.0g                      CD4 count


cd4_cat           byte      %11.0g       cd4_cat        CD4 count


hiv               byte      %12.0g       hiv            HIV status


hungry            byte      %12.0g       hungry         How often household members go hungry


lastalc           byte      %20.0g       lastalc        Last time used alcohol


n_adults          byte      %12.0g                      # adults in household


n_child           byte      %12.0g                      # people <18yo in household


n_support         byte      %12.0g                      # people supported financially


religion          byte      %10.0g       religion       Religion


sex               byte      %12.0g       sex            Sex


-------------
Sorted by: hiv

.
. /* Now get brief summary statistics for all the variables in the dataset. The
> full command we'll use is summarize, but we can shorten that to save typing. */
.
. sum

    Variable |        Obs        Mean    Std. Dev.       Min        Max
-------------+---------------------------------------------------------
         age |      3,387    31.72808    9.850006         18         80
      auditc |      3,244    .6390259    1.483416          0         12
         cd4 |        999    329.2332    266.1177          1       1932
     cd4_cat |        999    1.618619    .9377913          0          3
         hiv |      3,389    .2959575    .4565393          0          1
-------------+---------------------------------------------------------
      hungry |      3,321    3.603132    .7817434          1          4
     lastalc |      3,377    .9227125    .9101892          0          2
    n_adults |      3,384    1.664598    1.737786          0         12
     n_child |      3,382    1.850089    1.924868          0         12
   n_support |      3,372    3.907177    3.227549          0         20
-------------+---------------------------------------------------------
    religion |      3,387    3.781518    1.994997          1          8
         sex |      3,389    1.541753    .4983272          1          2

.
. /* Note how many people had CD4 counts.  Why is this less than the number of observations?
>
> To get more information about continuous variables, add the detail option to
> the summarize command. */
.
. sum age cd4, detail

                             Age
-------------------------------------------------------------
      Percentiles      Smallest
 1%           18             18
 5%           20             18
10%           21             18       Obs               3,387
25%           24             18       Sum of Wgt.       3,387

50%           30                      Mean           31.72808
                        Largest       Std. Dev.      9.850006
75%           38             75
90%           46             75       Variance       97.02261
95%           50             78       Skewness       1.030799
99%           60             80       Kurtosis       3.975972

                          CD4 count
-------------------------------------------------------------
      Percentiles      Smallest
 1%            5              1
 5%           14              2
10%           36              2       Obs                 999
25%          130              2       Sum of Wgt.         999

50%          283                      Mean           329.2332
                        Largest       Std. Dev.      266.1177
75%          463           1461
90%          659           1601       Variance       70818.64
95%          866           1804       Skewness       1.444705
99%         1182           1932       Kurtosis       6.518639

.
. /* The tabstat command also gives easy summaries. Here we use the stat() option
> to obtain the number of observations, mean, SD, minimum, 25th percentile,
> median, 75th percentile, and maximum. The format option controls rounding of
> the statistics. */
.
. tabstat age cd4, stat(n mean sd min p25 p50 p75 max) format(%8.3g)

   stats |       age       cd4
---------+--------------------
N   |      3387       999
mean   |      31.7       329
sd   |      9.85       266
min   |        18         1
p25   |        24       130
p50   |        30       283
p75   |        38       463
max   |        80      1932
---------+--------------------

.
. /* Adding the option col(statistics) transposes the results. Also, we can use
> the line break marker /// to tell stata that command continues after a hard
> return. This is handy if you want to use a relatively narrow do-file editor
> window. (Note that you can't use /// in the command window.)  */
.
. tabstat age cd4, stat(n mean sd min p25 p50 p75 max) ///
>         col(statistics) format(%8.3g)

variable   |         N      mean        sd       min       p25       p50       p75
-------------+----------------------------------------------------------------------
age   |      3387      31.7      9.85        18        24        30        38
cd4   |       999       329       266         1       130       283       463
-------------+----------------------------------------------------------------------

variable   |       max
-------------+----------
age   |        80
cd4   |      1932
-------------+----------

.
. /* Add the by() option to see the statistics stratified by HIV status. */
. tabstat age cd4, by(hiv) stat(n mean sd min p25 p50 p75 max) ///
>         col(statistics) format(%8.3g)

Summary for variables: age cd4
     by categories of: hiv (HIV status)

hiv   |         N      mean        sd       min       p25       p50       p75       max
-------+--------------------------------------------------------------------------------
HIV-   |      2384      31.3      10.5        18        23        28        38        80
       |         0         .         .         .         .         .         .         .
-------+--------------------------------------------------------------------------------
HIV+   |      1003      32.6      8.05        18        27        31        38        70
       |       999       329       266         1       130       283       463      1932
-------+--------------------------------------------------------------------------------
Total   |      3387      31.7      9.85        18        24        30        38        80
       |       999       329       266         1       130       283       463      1932
-------+--------------------------------------------------------------------------------

.
. /* The tabstat command can only be stratify by one by variable at a time; to
> stratify by >1 variables, use the table command.*/
. table hiv sex hungry, by(religion) c(mean age) format(%8.1f)

-----------+---------------------------------------------------------------------
Religion   |            How often household members go hungry and Sex
and HIV    | ---- Often ---    -- Sometimes -    --- Seldom ---    ---- Never ---
status     |   Male  Female      Male  Female      Male  Female      Male  Female
-----------+---------------------------------------------------------------------
Protestant |
      HIV- |   35.0    32.2      33.9    34.9      29.9    32.6      31.9    31.9
      HIV+ |   38.0    32.4      33.4    31.0      39.6    34.2      35.5    31.4
-----------+---------------------------------------------------------------------
Catholic   |
      HIV- |   33.5    38.7      29.2    34.2      32.5    33.4      31.3    31.2
      HIV+ |   33.7    36.3      33.9    29.8      33.8    31.1      34.1    30.5
-----------+---------------------------------------------------------------------
Saved      |
      HIV- |   24.5    34.4      30.1    28.4      26.1    32.4      29.5    30.3
      HIV+ |   38.0    27.0      36.0    30.6      31.7    32.4      35.8    30.9
-----------+---------------------------------------------------------------------
Moslem     |
      HIV- |   37.7    36.2      30.0    34.0      27.7    31.2      31.1    31.2
      HIV+ |   44.0    33.4      35.9    32.9      38.5    31.8      34.9    30.1
-----------+---------------------------------------------------------------------
Other      |
      HIV- |                     28.0    28.0      26.2    24.0      29.5    31.3
      HIV+ |           40.0      39.0              46.0    37.0      34.8    32.7
-----------+---------------------------------------------------------------------

.
. /* Exploratory data analysis using graphical means.  First make a histogram
> of the age distribution.  We can use the name option so that multiple graphs
> can be shown in the viewer.  The replace option is needed to prevent Stata
> from stopping if the graph already exists -- again, this is useful if you
> are rerunning a do-file to debug it. */
.
. histogram age, fcolor(blue) name(hist_age_00, replace)
(bin=35, start=18, width=1.7714286)

.
. /* Now add titles to the histogram. */
.
. histogram age, fcolor(blue) title(Histogram of age) ///
>         xtitle(Age) name(hist_age_01, replace)
(bin=35, start=18, width=1.7714286)

.
. /* By default, the Y-axis shows the "density", which sums to one for all bars.
> Instead we can get a histogram of frequencies. */
.
. histogram age, fcolor(blue) frequency ///
>         title(Histogram of age) xtitle(Age) name(hist_age_02, replace)
(bin=35, start=18, width=1.7714286)

.
. /* Now combine the two histograms into a single plot. Note that he histograms
> are identical in shape, and both show that age is somewhat skewed to the right. */
.
. graph combine hist_age_01 hist_age_02, name(hist_age_combined, replace)

.
. * save this as a PDF file
. graph export hist_age_combined.pdf, replace
(file /Users/eric/Work/MB/ATCR/200/Labs/hist_age_combined.pdf written in PDF format)

.
. /* Now we will change the number of bars in the histogram, leaving Stata to
> decide on their widths.
>
> We will do this using "loops." Specifically, the program will loop over the 4
> values of nb in the first line, pasting those values as so-called local
> variables `nb' into the histogram command in several places.  Learning to use
> loops can make your programming a lot less tedious. Recall that you will need
> to run the four lines of code defining the loop together. */
.
. foreach nb in 5 10 25 50 {
  2.         histogram age, fcolor(blue) percent  bin(`nb') ///
>                 title(`nb' bins) xtitle(Age) name(hanb`nb', replace)
  3.         }
(bin=5, start=18, width=12.4)
(bin=10, start=18, width=6.2)
(bin=25, start=18, width=2.48)
(bin=50, start=18, width=1.24)

. graph combine hanb5 hanb10 hanb25 hanb50, name(hist_age_bins, replace)

. graph export hist_age_bins.pdf, replace
(file /Users/eric/Work/MB/ATCR/200/Labs/hist_age_bins.pdf written in PDF format)

.
. /* As an alternative, we can vary the width of the bins, again using a loop. */
.
. foreach bw in 1 2 5 10 {
  2.         histogram age, fcolor(blue) percent width(`bw') ///
>                 title(bin width `bw') xtitle(Age) name(habw`bw', replace)
  3.         }
(bin=62, start=18, width=1)
(bin=31, start=18, width=2)
(bin=13, start=18, width=5)
(bin=7, start=18, width=10)

. graph combine habw1 habw2 habw5 habw10, name(hist_age_widths, replace)

. graph export hist_age_widths.pdf, replace
(file /Users/eric/Work/MB/ATCR/200/Labs/hist_age_widths.pdf written in PDF format)

.
. /* Box plots are also useful for assessing the distribution of continuous variables.  The box runs from the 25th to the 75th percentile, and includes
> a line at the median.  The so-called whiskers extend 1.5 IQRs from the bottom
> and top of the box, or to the most extreme value, whichever is closer to the
> box. Any outliers more than 1.5 IQRs from the box are plotted separately. */
.
. graph box age, name(box_age, replace)

.
. /* The boxplot also shows that age is somewhat right skewed, as we saw in the histograms.
>
> We can stratify the boxplot by up to three other factors. Note that the value
> labels we added to sex and hiv make this plot easier to understand */
.
. graph box age, over(sex) over(hiv) name(box_age_stratified, replace)

. graph export box_age_stratified.pdf, replace
(file /Users/eric/Work/MB/ATCR/200/Labs/box_age_stratified.pdf written in PDF format)

.
. /* Now we will use scatter plots to look at the association of a continuous
> outcome with a continuous predictor. */
.
. scatter n_support age, name(scatter_age, replace)

.
. /* Now we add a Lowess smooth through the data. Essentially this is a
> regression of n_support on age, without making the assumption that the
> regression line is linear. The bwidth option controls the smoothness of
> the line.  The default is 0.8, which means that Stata uses up to 80% of
> the data to determine the vertical location of the smooth at each value
> of the independent variable. The lineopts() option allows us to control
> the thickness of the line, and could also be used to reset its color and
> whether the line is solid, dashed, etc.
>
> We will again using loops to explore the bandwidth. Here, we have to define a
> second local variable for the bandwidth, based on looping variable pct,
> because the plot name can't include decimals. Note that the first line uses
> a somewhat different syntax than the previous loops.
> */
.
. foreach pct of numlist 80(-20)20 10 1 {
  2.         local bw = `pct'/100
  3.         lowess n_support age, ///
>                 bwidth(`bw') msize(tiny) lineopts(lwidth(thick)) title("") ///
>                 name(bw`pct', replace)
  4.         }

. graph combine bw80 bw60 bw40 bw20 bw10 bw1, name(lowess_combined, replace)

. graph export lowess_combined.pdf, replace
(file /Users/eric/Work/MB/ATCR/200/Labs/lowess_combined.pdf written in PDF format)

.
. /* Reducing the bandwith imposes less smoothness, but clearly the smallest  bandwidth is going too far!
>
> We can also focus on a certain age range by adding an if statement. */
.
. lowess n_support age if age<60, ///
>         msize(tiny) lineopts(lwidth(thick)) title("") ///
>         name(lowess_restricted, replace)

. graph export lowess_restricted.pdf, replace
(file /Users/eric/Work/MB/ATCR/200/Labs/lowess_restricted.pdf written in PDF format)

.
. /* Finally we can also use tabstat and box plots to look at the distribution
> of a continuous outcome when the independent variable is categorical. */
.
. tabstat n_support, by(hungry) stat(n mean sd min p25 p50 p75 max) format(%8.3g)

Summary for variables: n_support
     by categories of: hungry (How often household members go hungry)

hungry   |         N      mean        sd       min       p25       p50       p75
----------+----------------------------------------------------------------------
Often   |        77      4.38      3.36         0         1         4         6
Sometimes   |       378      4.18      3.18         0         2         4         6
Seldom   |       319       4.2      3.34         0         2         4         6
Never   |      2530      3.81      3.21         0         1         3         5
----------+----------------------------------------------------------------------
Total   |      3304      3.91      3.22         0         1         3         5
----------+----------------------------------------------------------------------

hungry   |       max
----------+----------
Often   |        13
Sometimes   |        16
Seldom   |        20
Never   |        20
----------+----------
Total   |        20
----------+----------

.
. graph box n_support, over(hungry) name(ins_alc, replace)

. graph export ins_alc.pdf, replace
(file /Users/eric/Work/MB/ATCR/200/Labs/ins_alc.pdf written in PDF format)

.
. /* Now close the log file.*/
.
. log close
      name:  <unnamed>
       log:  /Users/eric/Work/MB/ATCR/200/Labs/Lab 1 exercises 2018.smcl
  log type:  smcl
 closed on:   5 Jul 2018, 05:33:30
-------------