-------------
      name:  <unnamed>
       log:  /Users/eric/Work/MB/ATCR/200/Labs/lab 1.smcl
  log type:  smcl
 opened on:  18 Sep 2018, 06:54:15

.
. /* Read in the 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 Kampala, Uganda. */
.
. use vct.dta, clear

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

.
. /* First see what variables are in the dataset.  Check how many observations
> and variables the dataset includes. */
.
. describe

Contains data from vct.dta
  obs:         3,389
 vars:            13                          23 Jul 2018 10:11
 size:        71,169
-------------
              storage   display    value
variable name   type    format     label      variable label
-------------


age               double    %12.0g


auditc            byte      %12.0g


cd4               int       %11.0g


cd4_cat           byte      %11.0g


female            byte      %9.0g


hiv               byte      %12.0g


hungry            byte      %12.0g


lastalc           byte      %20.0g


n_adults          byte      %12.0g


n_child           byte      %12.0g


n_support         byte      %12.0g


religion          byte      %10.0g


sex               byte      %12.0g


-------------
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.73267    9.846937      17.52      79.85
      auditc |      3,244    .6390259    1.483416          0         12
         cd4 |        999    329.2332    266.1177          1       1932
     cd4_cat |        999    1.618619    .9377913          0          3
      female |      3,389    .5417527    .4983272          0          1
-------------+---------------------------------------------------------
         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

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

. label var auditc "Audit score, last 3 months"

. 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 female "female 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). Other easy-to-remember conventions will also work (e.g.,
> add an underscore to the variable name).*/
.
. label define yesno 0 "no" 1 "yes"

. label values female yesno

. label define sex 1 "men" 2 "women"

. label values sex sex

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

. label values hiv hiv

. label define cd4_cat 0 "<50" 1 "51-250" 2 "251-500" 3 ">500"

. label values cd4_cat cd4_cat

. label define hungry 1 "often" 2 "sometimes" 3 "seldom" 4 "never"

. label values hungry hungry

. label define lastalc 0 "never" 1 ">1 year ago" 2 "within last year"

. label values lastalc lastalc

. label define religion 1 "Protestant" 2 "Catholic" 3 "Saved" 4 "Moslem" 8 "Other"

. label values religion religion

.
. /* You can use these commands to get information about the values labels for
> a variable:*/
.
. codebook religion

-------------
religionReligion
-------------

                  type:  numeric (byte)
label  :  religion, but 2 nonmissing values are not labeled

                 range:  [1,8]                        units:  1
         unique values:  5                        missing .:  2/3,389

            tabulation:  Freq.   Numeric  Label
                         1,002         1  Protestant
                         1,047         4  Moslem
                           576         5
                           679         6
                            83         8  Other
                             2         .

. label list religion
religion:
           1 Protestant
           2 Catholic
           3 Saved
           4 Moslem
           8 Other

.
. /* Note the use of the indicator female for sex, coded 0 = no and 1 = yes.
> Indicators coded like this have several advantages:
>         * the mean value of the indicator gives the proportion of women in the dataset
>         * this is the coding required for binary outcomes in logistic models
>         * unlike a variable for sex, coded 1 = male 2 = female, the indicator does
>           not require us to remember the coding and value labels
>
> Now take another look at what we have. The value labels should show up in the
> browser. */
.
. des

Contains data from vct.dta
  obs:         3,389
 vars:            13                          23 Jul 2018 10:11
 size:        71,169
-------------
              storage   display    value
variable name   type    format     label      variable label
-------------


age               double    %12.0g                      Age


auditc            byte      %12.0g                      Audit score, last 3 months


cd4               int       %11.0g                      CD4 count


cd4_cat           byte      %11.0g       cd4_cat        CD4 count


female            byte      %9.0g        yesno          female sex


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


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

. sum

    Variable |        Obs        Mean    Std. Dev.       Min        Max
-------------+---------------------------------------------------------
         age |      3,387    31.73267    9.846937      17.52      79.85
      auditc |      3,244    .6390259    1.483416          0         12
         cd4 |        999    329.2332    266.1177          1       1932
     cd4_cat |        999    1.618619    .9377913          0          3
      female |      3,389    .5417527    .4983272          0          1
-------------+---------------------------------------------------------
         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

. browse

.
. /* Note that CD4 count is present only for the 30% of participants who are HIV+.
>
> ********************************************************************************
> Creating and modifying new variables
>
> To create a new variable, use the generate command; to modify it, use the
> replace command, and to drop a variable, use the drop command. In making
> variables, you will often need to use so-called logical operators:
>         equal:                          ==
>         not equal:              ~= or ~=
>         and:                                    &
>         or:                                             |
>         greater:                        > or >=
>         less:                                   < or <=
> Missing numeric values are denoted . (for missing character strings use "").
> Note that missing numeric values are stored as a very large number, (always
> the largest in the dataset).  As a result, if you are missing age values for
> 2 persons in the data set, and you create a new variable that includes all
> participants with age>30, you will include the 2 with missing data.
>
> Example: AUDIT-C is a scale measuring alcohol consumption, ranging from 0-12.
> For those not currently drinking (on this survey it is the past 3 months) the
> score is 0, for those drinking the score is 1-12.  To make a new variable that
> represents any alcohol consumption we can use this following code:  */
.
. gen alc_any=0 if auditc==0
(966 missing values generated)

. replace alc_any=1 if auditc>0 & auditc !=.
(821 real changes made)

. label var alc_any "Any alcohol use"

. label values alc_any yesno // value label yesno is already part of the dataset

. * check to make sure this worked
. tab auditc alc_any, missing

     Audit |
    score, |
    last 3 |         Any alcohol use
    months |        no        yes          . |     Total
-----------+---------------------------------+----------
         0 |     2,423          0          0 |     2,423
         1 |         0        363          0 |       363
         2 |         0        134          0 |       134
         3 |         0        118          0 |       118
         4 |         0        108          0 |       108
         5 |         0         44          0 |        44
         6 |         0         13          0 |        13
         7 |         0         10          0 |        10
         8 |         0          7          0 |         7
         9 |         0         15          0 |        15
        10 |         0          4          0 |         4
        11 |         0          3          0 |         3
        12 |         0          2          0 |         2
         . |         0          0        145 |       145
-----------+---------------------------------+----------
     Total |     2,423        821        145 |     3,389

.
. /* The second condition ensures that missings are kept as missing.  A better
> way to do this is to use the recode command, which will handle the missing
> values correctly: */
.
. drop alc_any

. recode auditc (0=0 "no") (2/max=1 "yes"), gen(alc_any)
(458 differences between auditc and alc_any)

. label var alc_any "Any alcohol use"

. * check to make sure this worked
. tab auditc alc_any, missing

     Audit |
    score, |
    last 3 |         Any alcohol use
    months |        no        yes          . |     Total
-----------+---------------------------------+----------
         0 |     2,423          0          0 |     2,423
         1 |         0        363          0 |       363
         2 |         0        134          0 |       134
         3 |         0        118          0 |       118
         4 |         0        108          0 |       108
         5 |         0         44          0 |        44
         6 |         0         13          0 |        13
         7 |         0         10          0 |        10
         8 |         0          7          0 |         7
         9 |         0         15          0 |        15
        10 |         0          4          0 |         4
        11 |         0          3          0 |         3
        12 |         0          2          0 |         2
         . |         0          0        145 |       145
-----------+---------------------------------+----------
     Total |     2,423        821        145 |     3,389

.
. /* We can also use recode and replace to make a more complex indicator for
> unhealthy alcohol consumption, using different cutoff for men (>=4) and
> women (>=3): */
.
. recode auditc (min/2=0 "no") (3/max=1 "yes"), gen(alc_unh)
(821 differences between auditc and alc_unh)

. replace alc_unh = 0 if female==0 & auditc==3
(69 real changes made)

. label var alc_unh "Unhealthy alcohol use"

. * check to make sure this worked
. bysort female: tab auditc alc_unh, m

-------------
-> female = no

     Audit |
    score, |
    last 3 |      Unhealthy alcohol use
    months |        no        yes          . |     Total
-----------+---------------------------------+----------
         0 |     1,042          0          0 |     1,042
         1 |       166          0          0 |       166
         2 |        57          0          0 |        57
         3 |        69          0          0 |        69
         4 |         0         71          0 |        71
         5 |         0         34          0 |        34
         6 |         0          9          0 |         9
         7 |         0          8          0 |         8
         8 |         0          5          0 |         5
         9 |         0         12          0 |        12
        10 |         0          2          0 |         2
        11 |         0          2          0 |         2
        12 |         0          2          0 |         2
         . |         0          0         74 |        74
-----------+---------------------------------+----------
     Total |     1,334        145         74 |     1,553

-------------
-> female = yes

     Audit |
    score, |
    last 3 |      Unhealthy alcohol use
    months |        no        yes          . |     Total
-----------+---------------------------------+----------
         0 |     1,381          0          0 |     1,381
         1 |       197          0          0 |       197
         2 |        77          0          0 |        77
         3 |         0         49          0 |        49
         4 |         0         37          0 |        37
         5 |         0         10          0 |        10
         6 |         0          4          0 |         4
         7 |         0          2          0 |         2
         8 |         0          2          0 |         2
         9 |         0          3          0 |         3
        10 |         0          2          0 |         2
        11 |         0          1          0 |         1
         . |         0          0         71 |        71
-----------+---------------------------------+----------
     Total |     1,655        110         71 |     1,836


.
. /* Note the abbreviated name for the unhealthy alcohol use indicator.  Stata
> does allow long, explicit variable names, with some limits, but they can be
> tedious to use in programming; better to encode the extra information in the
> variable label.   Note also that all the alcohol variables we have made start
> with the same stem alc_, also with a view towards making them easier to remember.
>
> Finally, we can use recode to categorize the continuous variable, age.  In
> doing this, we can use the ordering of the elements of the command to control
> whether the values at each cutpoint are put into the lower or upper category.
> Specifically, values equal to a cutpoint get put into the class defined by the
> element where the cutpoint is first encountered. We will first code age
> categories <=20, >20 to <=30, >30 to <=40, >40 to <=50, and >50. */
.
. qui recode age (min/20=1 "<=20") (20/30=2 ">20 to 30") (30/40=3 ">30 to 40") ///
>         (40/50=4 ">40 to 50") (50/max=5 ">50"), gen(age_cat)

. label var age_cat "Age"

. * check this
. tabstat age, by(age_cat) s(min max)

Summary for variables: age
     by categories of: age_cat (Age)

age_cat   |       min       max
----------+--------------------
<=20   |     17.52        20
>20 to 30   |     20.01        30
>30 to 40   |     30.01        40
>40 to 50   |     40.01     49.99
>50   |     50.01     79.85
----------+--------------------
Total   |     17.52     79.85
----------+--------------------

.
. /* To include the cutpoints in the upper category, we can reverse the order of
> the elements of the command: */
.
. drop age_cat

. qui recode age  (50/max=5 "50+") (40/50=4 "40 to <50") (30/40=3 "30 to <40") ///
>         (20/30=2 "20 to <30") (min/20=1 "<20") , gen(age_cat)

. label var age_cat "Age"

. * check this
. tabstat age, by(age_cat) s(min max)

Summary for variables: age
     by categories of: age_cat (Age)

age_cat   |       min       max
----------+--------------------
<20   |     17.52     19.98
20 to <30   |        20     29.99
30 to <40   |        30     39.99
40 to <50   |        40     49.99
50+   |     50.01     79.85
----------+--------------------
Total   |     17.52     79.85
----------+--------------------

.
. /* Put the variables in alphabetical order, compress anything stored wastefully,
> get a final look at the dataet, and save the modified dataset. */
.
. compress
  variable age_cat was double now byte
  (23,723 bytes saved)

. aorder

.
. /* Finally, save the revised dataset.  It's good practice never to write over
> the input dataset. */
.
. save vct_revised, replace
file vct_revised.dta saved

.
. /* Now we will use the revised dataset to do more exploration. */
.
. use vct_revised.dta, clear

.
. /* Briefly remind yourself of what is in the dataset. Check how many observations
> and variables the dataset includes. */
.
. describe

Contains data from vct_revised.dta
  obs:         3,389
 vars:            16                          18 Sep 2018 06:54
 size:        81,336
-------------
              storage   display    value
variable name   type    format     label      variable label
-------------


age               double    %12.0g                      Age


age_cat           byte      %9.0g        age_cat        Age


alc_any           byte      %9.0g        alc_any        Any alcohol use


alc_unh           byte      %9.0g        alc_unh        Unhealthy alcohol use


auditc            byte      %12.0g                      Audit score, last 3 months


cd4               int       %11.0g                      CD4 count


cd4_cat           byte      %11.0g       cd4_cat        CD4 count


female            byte      %9.0g        yesno          female sex


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


-------------
Sorted by: female

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

    Variable |        Obs        Mean    Std. Dev.       Min        Max
-------------+---------------------------------------------------------
         age |      3,387    31.73267    9.846937      17.52      79.85
     age_cat |      3,387    2.653085    .9790142          1          5
     alc_any |      3,244    .2530826    .4348449          0          1
     alc_unh |      3,244    .0786067    .2691653          0          1
      auditc |      3,244    .6390259    1.483416          0         12
-------------+---------------------------------------------------------
         cd4 |        999    329.2332    266.1177          1       1932
     cd4_cat |        999    1.618619    .9377913          0          3
      female |      3,389    .5417527    .4983272          0          1
         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.04          17.52
 5%        19.65          17.56
10%        21.15          17.57       Obs               3,387
25%        24.25           17.6       Sum of Wgt.       3,387

50%        29.53                      Mean           31.73267
                        Largest       Std. Dev.      9.846937
75%        37.73          74.62
90%        45.89          75.33       Variance       96.96217
95%        50.42          78.29       Skewness       1.029901
99%        60.01          79.85       Kurtosis       3.974463

                          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

.
.
. ********************************************************************************
.
. /* Exploratory analysis using descriptive statistics
>
> The tabstat command 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   |      17.5         1
p25   |      24.3       130
p50   |      29.5       283
p75   |      37.7       463
max   |      79.9      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      17.5      24.3      29.5      37.7
cd4   |       999       329       266         1       130       283       463
-------------+----------------------------------------------------------------------

variable   |       max
-------------+----------
age   |      79.9
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      17.5      23.4      28.3      37.5      79.9
       |         0         .         .         .         .         .         .         .
-------+--------------------------------------------------------------------------------
HIV+   |      1003      32.7      8.05      17.6      26.8      31.4      38.1        70
       |       999       329       266         1       130       283       463      1932
-------+--------------------------------------------------------------------------------
Total   |      3387      31.7      9.85      17.5      24.3      29.5      37.7      79.9
       |       999       329       266         1       130       283       463      1932
-------+--------------------------------------------------------------------------------

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

-----------+-----------------------------------------------------
Religion   | How often household members go hungry and female sex
and HIV    | -- often -    - someti -    - seldom -    -- never -
status     |   no   yes      no   yes      no   yes      no   yes
-----------+-----------------------------------------------------
Protestant |
      HIV- | 34.9  32.3    33.9  34.9    29.8  32.5    31.9  31.9
      HIV+ | 38.3  32.6    33.5  31.1    39.5  34.2    35.4  31.4
-----------+-----------------------------------------------------
Moslem     |
      HIV- | 33.4  38.7    29.2  34.2    32.5  33.3    31.3  31.2
      HIV+ | 34.0  36.2    33.9  29.9    33.8  31.1    34.1  30.5
-----------+-----------------------------------------------------
5          |
      HIV- | 24.9  34.4    30.1  28.4    26.1  32.5    29.5  30.3
      HIV+ | 38.3  26.9    35.9  30.5    31.5  32.4    35.8  30.9
-----------+-----------------------------------------------------
6          |
      HIV- | 37.8  36.2    30.0  34.0    27.7  31.3    31.1  31.2
      HIV+ | 43.6  33.4    35.9  33.1    38.7  31.8    34.9  30.0
-----------+-----------------------------------------------------
Other      |
      HIV- |               27.9  28.0    26.1  24.2    29.5  31.3
      HIV+ |       39.6    39.5          46.5  37.1    35.0  32.8
-----------+-----------------------------------------------------

.
. * To see how categorical variables relate to each other, use the tab command.
. tab age_cat hiv, row

+----------------+
| Key            |
+----------------+
|   frequency    |
| row percentage |
+----------------+

           |      HIV status
       Age |      HIV-       HIV+ |     Total
-----------+----------------------+----------
       <20 |       196         24 |       220
           |     89.09      10.91 |    100.00
-----------+----------------------+----------
 20 to <30 |     1,150        412 |     1,562
           |     73.62      26.38 |    100.00
-----------+----------------------+----------
 30 to <40 |       576        394 |       970
           |     59.38      40.62 |    100.00
-----------+----------------------+----------
 40 to <50 |       302        141 |       443
           |     68.17      31.83 |    100.00
-----------+----------------------+----------
       50+ |       160         32 |       192
           |     83.33      16.67 |    100.00
-----------+----------------------+----------
     Total |     2,384      1,003 |     3,387
           |     70.39      29.61 |    100.00

. tab female hiv, row

+----------------+
| Key            |
+----------------+
|   frequency    |
| row percentage |
+----------------+

           |      HIV status
female sex |      HIV-       HIV+ |     Total
-----------+----------------------+----------
        no |     1,175        378 |     1,553
           |     75.66      24.34 |    100.00
-----------+----------------------+----------
       yes |     1,211        625 |     1,836
           |     65.96      34.04 |    100.00
-----------+----------------------+----------
     Total |     2,386      1,003 |     3,389
           |     70.40      29.60 |    100.00

. tab hungry hiv, row

+----------------+
| Key            |
+----------------+
|   frequency    |
| row percentage |
+----------------+

 How often |
 household |
members go |      HIV status
    hungry |      HIV-       HIV+ |     Total
-----------+----------------------+----------
     often |        49         30 |        79
           |     62.03      37.97 |    100.00
-----------+----------------------+----------
 sometimes |       246        134 |       380
           |     64.74      35.26 |    100.00
-----------+----------------------+----------
    seldom |       209        112 |       321
           |     65.11      34.89 |    100.00
-----------+----------------------+----------
     never |     1,841        700 |     2,541
           |     72.45      27.55 |    100.00
-----------+----------------------+----------
     Total |     2,345        976 |     3,321
           |     70.61      29.39 |    100.00

. tab alc_any hiv, row

+----------------+
| Key            |
+----------------+
|   frequency    |
| row percentage |
+----------------+

       Any |
   alcohol |      HIV status
       use |      HIV-       HIV+ |     Total
-----------+----------------------+----------
        no |     1,753        670 |     2,423
           |     72.35      27.65 |    100.00
-----------+----------------------+----------
       yes |       540        281 |       821
           |     65.77      34.23 |    100.00
-----------+----------------------+----------
     Total |     2,293        951 |     3,244
           |     70.68      29.32 |    100.00

.
. /* Taking advantage of the 0-1 coding of the HIV indicator, we can use the
> table command get HIV prevalence by several other variables at once; here we
> do just age and sex. */
.
. table age_cat female, c(mean hiv) format(%8.2f)

----------+-----------
          | female sex
      Age |   no   yes
----------+-----------
      <20 | 0.03  0.18
20 to <30 | 0.15  0.35
30 to <40 | 0.40  0.42
40 to <50 | 0.35  0.29
      50+ | 0.20  0.13
----------+-----------

.
. ********************************************************************************
.
. /* We can also use the Stata means command to calculate confidence intervals
> for the means of continuous variables, both overall and within strata. */
.
. mean age, cformat(%8.2f)

Mean estimation                   Number of obs   =      3,387

-------------+------------------------------------------------
             |       Mean    Std. Err.      [95% Con f. Interval]
-------------+------------------------------------------------
         age |       31.73        0.17          31.40        32.06
-------------+------------------------------------------------

. mean age, over(female hiv) cformat(%8.2f)

Mean estimation                   Number of obs   =      3,387

         Over: female hiv
    _subpop_1: no HIV-
    _subpop_2: no HIV+
    _subpop_3: yes HIV-
    _subpop_4: yes HIV+

-------------+------------------------------------------------
         Over |       Mean    Std. Err.      [95% Con f. Interval]
-------------+------------------------------------------------
age          |
   _subpop_1 |       30.99        0.31          30.39        31.60
   _subpop_2 |       35.17        0.39          34.39        35.94
   _subpop_3 |       31.68        0.30          31.10        32.27
   _subpop_4 |       31.15        0.32          30.53        31.77
-------------+------------------------------------------------

.
. * To get, say, 90% confidence intervals, we can use the level option
. mean age, over(female hiv) level(90) cformat(%8.2f)

Mean estimation                   Number of obs   =      3,387

         Over: female hiv
    _subpop_1: no HIV-
    _subpop_2: no HIV+
    _subpop_3: yes HIV-
    _subpop_4: yes HIV+

-------------+------------------------------------------------
         Over |       Mean    Std. Err.      [90% Con f. Interval]
-------------+------------------------------------------------
age          |
   _subpop_1 |       30.99        0.31          30.48        31.50
   _subpop_2 |       35.17        0.39          34.52        35.81
   _subpop_3 |       31.68        0.30          31.19        32.18
   _subpop_4 |       31.15        0.32          30.63        31.67
-------------+------------------------------------------------

.
. /* The means command could also be used to get means and confidence intervals
> for binary variables coded 0-1, but a better approach is the ci prop command,
> with the exact option.  This avoids confidence limits less than 0 or more
> than 1. */
.
. ci prop hiv, exact

                                                         -- Binomial Exact --
    Variable |        Obs  Proportion    Std. Err.       [95% Conf. Interval]
-------------+---------------------------------------------------------------
         hiv |      3,389    .2959575    .0078411        .2806231    .3116447

. ci prop hiv, exact level(90)

                                                         -- Binomial Exact --
    Variable |        Obs  Proportion    Std. Err.       [90% Conf. Interval]
-------------+---------------------------------------------------------------
         hiv |      3,389    .2959575    .0078411        .2830462    .3091298

.
. /* The ci command does not allow stratification, but we can get around that
> using the bysort prefix: */
.
. bysort female: ci prop hiv, exact

-------------
-> female = no

                                                         -- Binomial Exact --
    Variable |        Obs  Proportion    Std. Err.       [95% Conf. Interval]
-------------+---------------------------------------------------------------
         hiv |      1,553    .2433999    .0108895        .2222362    .2655381

-------------
-> female = yes

                                                         -- Binomial Exact --
    Variable |        Obs  Proportion    Std. Err.       [95% Conf. Interval]
-------------+---------------------------------------------------------------
         hiv |      1,836    .3404139    .0110587        .3187364    .3626026

.
. * The ci command also works for continous and Poisson-distributed count varibles
. ci means age

    Variable |        Obs        Mean    Std. Err.       [95% Conf. Interval]
-------------+---------------------------------------------------------------
         age |      3,387    31.73267    .1691974        31.40093    32.06441

. ci means n_adults, poisson

                                                         -- Poisson  Exact --
    Variable |   Exposure        Mean    Std. Err.       [95% Conf. Interval]
-------------+---------------------------------------------------------------
    n_adults |       3384    1.664598    .0221789        1.621409    1.708646

.
. ********************************************************************************
.
. /* 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=17.52, width=1.7808571)

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

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

.
. /* 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)

.
. /* To save this as a PDF file, we can use the following command.  Other options
> include saving as png, wmf, tif, eps (encapsulated postscript), and svg
> (scalable vector graphics), which some journals prefer.  See help graph export
> for more information. */
.
. 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 nbins (i.e., # of bins) in the first line, pasting those values as
> so-called local variables `nbins' into the histogram command in several places.
> Learning to use loops can make programming a lot less tedious. You will need
> to run the four lines of code defining the loop together. */
.
. foreach nbins in 5 10 25 50 {
  2.         histogram age, fcolor(blue) percent  bin(`nbins') ///
>                 title("`nbins' bins") xtitle(Age) name(hanb`nbins', replace)
  3.         }
(bin=5, start=17.52, width=12.466)
(bin=10, start=17.52, width=6.233)
(bin=25, start=17.52, width=2.4932)
(bin=50, start=17.52, width=1.2466)

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

.
. /* 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=63, start=17.52, width=1)
(bin=32, start=17.52, width=2)
(bin=13, start=17.52, width=5)
(bin=7, start=17.52, width=10)

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

.
. ********************************************************************************
.
. /* 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 in lab 1 to sex and hiv make this plot easier to understand */
.
. graph box age, over(sex) over(hiv) name(box_age_stratified, replace)

.
. /* Finally, Q-Q (quantile-quantile) plots are a good way to assess Normality.
> We will use a loop to look at them within groups, using some new commands to
> control the loops and title the plots: */
.
. qui levelsof sex, local(sl) // extracts distinct values of sex

. foreach i in `sl' {
  2.         local sex : label(sex) `i' // gets value label for level i of sex
  3.         forvalues j = 0/1 {
  4.                 local stat : label(hiv) `j'  // gets value label for level j of hiv
  5.                 qnorm age if sex==`i' & hiv==`j', ///
>                         title("`stat' `sex'") name(qq_age`i'`j', replace)
  6.                 }
  7.         }

. graph combine qq_age10 qq_age11 qq_age20 qq_age21, name(qq_age, replace)

.
. /* With the exception of HIV positive men, these plots all show the same right
> skewness we observed in the box plots and histograms. */
.
. ********************************************************************************
.
. /* 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 the nearest 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. The looping variable is a
> percentage, because the plot name can't include decimals. This means that we
> have to use some arcane code to convert `pct' to a proportion for use in the
> bw() option in the lowess command.  Note that the first line of the loop uses
> a somewhat different syntax than the previous loops.  */
.
. foreach pct of numlist 80(-20)20 10 1 {
  2.         lowess n_support age, ///
>                 bwidth(`=`pct'/100') msize(tiny) lineopts(lwidth(medthick)) title("") ///
>                 name(bw`pct', replace)
  3.         }

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

.
. /* 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(medthick)) title("") ///
>         name(lowess_restricted, replace)

.
. /* Finally, we can use tabstat and box plots to look at the conditional
> 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)

.
. /* Now close the log file.*/
.
. log close
      name:  <unnamed>
       log:  /Users/eric/Work/MB/ATCR/200/Labs/lab 1.smcl
  log type:  smcl
 closed on:  18 Sep 2018, 06:54:49
-------------