-------------
      name:  <unnamed>
       log:  /Users/eric/Work/MB/ATCR/200/Labs/lab 4.smcl
  log type:  smcl
 opened on:  20 Sep 2018, 11:29:51

.
. /* In this lab, we first show how to use various "immediate" commands in Stata,
> which are useful for running hypothesis tests using summary data (e.g., sample
> sizes, means, standard deviations, counts or proportions).  Then we will show
> how to run these tests using baseline data from the 2763 participants in the
> HERS RCT of hormone therapy as secondary prevention for coronary heart disease.
>
> ********************************************************************************
> Tests using immediate commands
>
> For the population of middle-aged men who later develop diabetes, baseline BMI
> has unknown mean and standard deviation(SD). In a sample of 58 men selected
> from this population, mean BMI is 25.0 kg/m2, with SD 2.7.
>
> First, calculate a 95% confidence interval for the mean. Begin by checking the
> help for the cii command. */
.
. help cii

.
. /* We have to specify the means option.  The arguments, in order, are the sample
> size, the mean, and the SD. */
.
. cii means 58 25 2.7

    Variable |        Obs        Mean    Std. Err.       [95% Conf. Interval]
-------------+---------------------------------------------------------------
             |         58          25    .3545274        24.29007    25.70993

.
. /* Now test the null hypothesis that the mean in this population is 24.0 kg/m2,
> the known mean of the population of men who do not develop diabetes.  The
> immediate commands are restricted to parametric tests, so we are implicitly
> assuming the BMI is reasonably normal, which is of course impossible to check
> from the summary data. */
.
. help ttesti

.
. /* We need the 1-sample t-test, comparing the observed sample mean to a known value.  The last argument for the command is the null value of the mean. */
.
. ttesti 58 25 2.7 24

One-sample t test
---------+--------------------------------------------------------------------
         |     Obs        Mean    Std. Err.   Std. Dev.   [95% Conf. Interval]
---------+--------------------------------------------------------------------
       x |      58          25    .3545274         2.7    24.29007    25.70993
---------+--------------------------------------------------------------------
    mean = mean(x)                                                t =   2.8207
Ho: mean = 24                                    degrees of freedom =       57

    Ha: mean < 24               Ha: mean != 24                 Ha: mean > 24
 Pr(T < t) = 0.9967         Pr(|T| > |t|) = 0.0066          Pr(T > t) = 0.0033

.
. /* So reasonably strong evidence that mean BMI among men who go on to develop
> diabetes is higher than the population mean among men who do not.
>
> **********
>
> In a study investigating the effect of seatbelt use on morbidity and mortality
> among pediatric victims of motor vehicles accidents, two random samples were
> selected, one of children wearing seat belts at the time of the accident
> (N=123), and the other of children who were not (N=290). In the first sample,
> 3 children died, while in the second, 13 did. Is there evidence that these
> two proportions are different?
>
> We can do this using the prtesti command. */
.
. help prtesti

.
. /* We need the immediate two-sample test of proportions with the count option. */
.
. prtesti 123 3 290 13, count

Two-sample test of proportions                     x: Number of obs =      123
                                                   y: Number of obs =      290
-------------+----------------------------------------------------------------
             |       Mean   Std. Err.      z    P>|z|     [95% Conf. Interval]
-------------+----------------------------------------------------------------
           x |   .0243902   .0139089                     -.0028707    .0516512
           y |   .0448276   .0121511                      .0210119    .0686432
-------------+----------------------------------------------------------------
        diff |  -.0204373   .0184691                     -.0566361    .0157614
             |  under Ho:   .0207648    -0.98   0.325
-------------+----------------------------------------------------------------
        diff = prop(x) - prop(y)                                  z =  -0.9842
    Ho: diff = 0

    Ha: diff < 0                 Ha: diff != 0                 Ha: diff > 0
 Pr(Z < z) = 0.1625         Pr(|Z| > |z|) = 0.3250          Pr(Z > z) = 0.8375

.
. /* We could also use csi command, which has the advantage of providing exact
> tests, which could be needed, since the number of deaths among children wearing
> seat belts was <5. It also provides confidence intervals for the risk difference
> and risk ratio.  In this case, the required input is the number of cases and
> non-cases, rather than the number of cases and sample size. We have set this
> up to treat seat belt use as exposure. */
.
. help csi

.
. csi 3 13 120 277, exact

                 |   Exposed   Unexposed  |      Total
-----------------+------------------------+------------
           Cases |         3          13  |         16
        Noncases |       120         277  |        397
-----------------+------------------------+------------
           Total |       123         290  |        413
                 |                        |
            Risk |  .0243902    .0448276  |   .0387409
                 |                        |
                 |      Point estimate    |    [95% Conf. Interval]
                 +------------------------+------------------------
 Risk difference |        -.0204373       |   -.0566361    .0157614
      Risk ratio |         .5440901       |    .1578389    1.875545
 Prev. frac. ex. |         .4559099       |   -.8755455    .8421611
 Prev. frac. pop |         .1357795       |
                 +------------------------+------------------------
                                  1-sided Fisher's exact P = 0.2468
                                  2-sided Fisher's exact P = 0.4126

.
. /* At least in this case, the chi-square and exact tests give qualitatively
> similar results. The difference may be NS, but the CIs for the risk difference
> and risk ratio are wide. A puzzing result -- are seatbelts ineffective for
> young children, or was the sample size too small?
>
> ***********
>
> Sixty-three adult males with coronary artery disease were stress-tested riding
> a stationary bicycle until the onset of chest pain. After resting, participants
> rode again until the repeat onset of angina, and the percent decrease in time
> to onset of pain was recorded. This procedure was carried out twice, first with
> clean air, and then, on a subsequent visit, when carbon monoxide had been mixed
> in the air. The sample size was 63 and the mean decreases in time to second
> onset were 3.35% and 9.63% at the first and second visits, and the standard
> deviation of the difference, 9.63-3.35=6.28%, was 20.3%. The research question
> is whether carbon monoxide affected the decrease in time to onset of chest
> pain on the second trial.
>
> We can test the difference using a paired t-test. */
.
. help ttesti

.
. /* We need the immediate for of the one sample test, with the null hypothesis
> (the last argument to the command) of zero. */
.
. ttesti 63 6.28 20.3 0

One-sample t test
---------+--------------------------------------------------------------------
         |     Obs        Mean    Std. Err.   Std. Dev.   [95% Conf. Interval]
---------+--------------------------------------------------------------------
       x |      63        6.28     2.55756        20.3    1.167511    11.39249
---------+--------------------------------------------------------------------
    mean = mean(x)                                                t =   2.4555
Ho: mean = 0                                     degrees of freedom =       62

    Ha: mean < 0                 Ha: mean != 0                 Ha: mean > 0
 Pr(T < t) = 0.9916         Pr(|T| > |t|) = 0.0169          Pr(T > t) = 0.0084

.
. /********************************************************************************
> Tests using actual data
>
> Now we will run analogous tests using actual rather than summarized data.
> Begin by reading in the HERS data. */
.
. use hers.dta, clear

.
. describe

Contains data from hers.dta
  obs:         2,763
 vars:            46                          23 Jul 2018 10:25
 size:       237,618
-------------
              storage   display    value
variable name   type    format     label      variable label
-------------


ptid              int       %9.0g                       anonymized participant identifier


acei              byte      %9.0g        noyes          use of ace-inhibitors


age               byte      %9.0g                       age (years)


anyllm            byte      %9.0g        noyes          use of any lipid-lowering agent


aspirin           byte      %9.0g        noyes          use of aspirin


betab             byte      %9.0g        noyes          use of beta-blockers


bmi               float     %9.0g                       bmi (kg/m^2)


ccb               byte      %9.0g        noyes          use of calcium channel blockers


crcl              float     %9.0g                       creatinine clearance (ml/min)


creat             float     %9.0g                       creatinine (mg/dl)


dbp               int       %9.0g                       diastolic blood pressure


diabetes          byte      %9.0g        noyes          diabetes


digitali          byte      %9.0g        noyes          use of digitalis


dmpills           byte      %9.0g        noyes          oral dm medication by self-report


drinkany          byte      %9.0g        noyes          any current alcohol consumption


drnkspwk          float     %9.0g                       average drinks per week


educyrs           byte      %9.0g                       years of education


exercise          byte      %9.0g        noyes          exercise at least 3 times per week


globrat           byte      %9.0g        globrat        self-reported health


glucose           int       %9.0g                       fasting glucose (mg/dl)


hdl0              int       %9.0g                       baseline hdl cholesterol (mg/dl)


hdl1              int       %9.0g                       year 1 hdl cholesterol (mg/dl)


ht                byte      %15.0g       ht             random assignment to hormone therapy


htnmeds           byte      %9.0g        noyes          anti-hypertensive use


insulin           byte      %9.0g        noyes          insulin use by self-report


ldl0              float     %9.0g                       baseline ldl cholesterol (mg/dl)


ldl1              float     %9.0g                       year 1 ldl cholesterol (mg/dl)


lpa0              float     %9.0g                       baseline Lp(a) (mg/dl)


lpa1              float     %9.0g                       year 1 Lp(a) (mg/dl)


medcond           byte      %9.0g                       other serious conditions by self-report


nonwhite          byte      %9.0g        noyes          nonwhite race/ethnicity


pcabg             byte      %9.0g        noyes          pre-randomization CABG


physact           byte      %20.0g       physact        comparative physical activity


pmi               byte      %9.0g        noyes          pre-randomization MI


poorfair          byte      %9.0g        noyes          poor/fair self-reported health


pptca             byte      %9.0g        noyes          pre-randomization PTCA


raceth            byte      %20.0g       raceth         race/ethnicity


sbp               int       %9.0g                       systolic blood pressure


smoking           byte      %9.0g        noyes          current smoker


statins           byte      %9.0g        noyes          use of statins


tchol0            int       %9.0g                       baseline total cholesterol (mg/dl)


tchol1            int       %9.0g                       year 1 total cholesterol (mg/dl)


tgl0              int       %9.0g                       baseline triglycerides (mg/dl)


tgl1              int       %9.0g                       year 1 triglycerides (mg/dl)


weight            float     %9.0g                       weight (kg)


whr               float     %9.0g                       waist/hip ratio


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

. sum

    Variable |        Obs        Mean    Std. Dev.       Min        Max
-------------+---------------------------------------------------------
        ptid |      2,763        1382    797.7537          1       2763
        acei |      2,763    .1758958    .3808004          0          1
         age |      2,763    66.64857    6.653081         44         79
      anyllm |      2,763    .4828085    .4997948          0          1
     aspirin |      2,763    .7915309    .4062874          0          1
-------------+---------------------------------------------------------
       betab |      2,763    .3260948    .4688673          0          1
         bmi |      2,758    28.57925    5.517783      15.21      54.13
         ccb |      2,763    .5479551    .4977851          0          1
        crcl |      2,760    61.30613    19.60505    6.18408   167.5208
       creat |      2,762    1.084432    .3163941         .4        9.3
-------------+---------------------------------------------------------
         dbp |      2,762     73.1517    9.705057         45        102
    diabetes |      2,763    .2645675    .4411825          0          1
    digitali |      2,763    .0998914    .2999095          0          1
     dmpills |      2,763    .0962722    .2950175          0          1
    drinkany |      2,761    .3915248    .4881797          0          1
-------------+---------------------------------------------------------
    drnkspwk |      2,763    1.355169    3.652089          0       49.5
     educyrs |      2,761    12.66932    2.645689          1         20
    exercise |      2,763    .3865364     .487044          0          1
     globrat |      2,760    3.063406     .842833          1          5
     glucose |      2,763    112.1516    36.84941         29        298
-------------+---------------------------------------------------------
        hdl0 |      2,752    50.26199    13.21566         14        130
        hdl1 |      2,608    51.78183    13.87166         14        124
          ht |      2,763    .4994571    .5000902          0          1
     htnmeds |      2,763    .8215708     .382943          0          1
     insulin |      2,763    .0988056    .2984549          0          1
-------------+---------------------------------------------------------
        ldl0 |      2,752    145.0385    37.80322       36.8      393.4
        ldl1 |      2,608    132.3915    39.02262        -20      450.2
        lpa0 |      2,759    33.71443    32.61081          0        236
        lpa1 |      2,612    30.95115    33.28073          0      201.9
     medcond |      2,763    .3720594    .4834416          0          1
-------------+---------------------------------------------------------
    nonwhite |      2,763    .1129207    .3165532          0          1
       pcabg |      2,763    .4096996    .4918673          0          1
     physact |      2,763    3.200145    1.083468          1          5
         pmi |      2,763    .5099529    .4999914          0          1
    poorfair |      2,760     .240942    .4277327          0          1
-------------+---------------------------------------------------------
       pptca |      2,763    .4303294    .4952118          0          1
      raceth |      2,763    1.146942     .439843          1          3
         sbp |      2,763    135.0695    19.02781         83        224
     smoking |      2,763    .1302932    .3366866          0          1
     statins |      2,763    .3731451    .4837277          0          1
-------------+---------------------------------------------------------
      tchol0 |      2,759    228.5799    41.01411        110        465
      tchol1 |      2,613    219.2346    41.10315         92        535
        tgl0 |      2,759    166.1493    63.51077         31        476
        tgl1 |      2,613    175.7976    84.45064         31       1016
      weight |      2,761     72.7305    14.70306       37.5      131.6
-------------+---------------------------------------------------------
         whr |      2,760     .869979    .0795729       .624      1.218

.
. /* First compare baseline LDL values in smokers and non-smokers using the
> t-test. */
.
. ttest ldl0, by(smoking)

Two-sample t test with equal variances
---------+--------------------------------------------------------------------
Variable |     Obs        Mean    Std. Err.   Std. Dev.   [95% Conf. Interval]
---------+--------------------------------------------------------------------
      no |   2,397    144.3985    .7639244     37.4011    142.9005    145.8965
     yes |     355      149.36    2.133767    40.20324    145.1635    153.5565
---------+--------------------------------------------------------------------
combined |   2,752    145.0385    .7206171    37.80322    143.6255    146.4515
---------+--------------------------------------------------------------------
    diff |           -4.961502    2.148142               -9.173636   -.7493675
---------+--------------------------------------------------------------------
    diff = mean(no) - mean(yes)                                   t =  -2.3097
Ho: diff = 0                                     degrees of freedom =     2750

    Ha: diff < 0                 Ha: diff != 0                 Ha: diff > 0
 Pr(T < t) = 0.0105         Pr(|T| > |t|) = 0.0210          Pr(T > t) = 0.9895

.
. /* We can check whether the results are sensitive to allowing for unequal
> variance in smoker and non-smokers. */
.
. ttest ldl0, by(smoking) unequal

Two-sample t test with unequal variances
---------+--------------------------------------------------------------------
Variable |     Obs        Mean    Std. Err.   Std. Dev.   [95% Conf. Interval]
---------+--------------------------------------------------------------------
      no |   2,397    144.3985    .7639244     37.4011    142.9005    145.8965
     yes |     355      149.36    2.133767    40.20324    145.1635    153.5565
---------+--------------------------------------------------------------------
combined |   2,752    145.0385    .7206171    37.80322    143.6255    146.4515
---------+--------------------------------------------------------------------
    diff |           -4.961502    2.266394               -9.415545   -.5074582
---------+--------------------------------------------------------------------
    diff = mean(no) - mean(yes)                                   t =  -2.1892
Ho: diff = 0                     Satterthwaite's degrees of freedom =  449.474

    Ha: diff < 0                 Ha: diff != 0                 Ha: diff > 0
 Pr(T < t) = 0.0145         Pr(|T| > |t|) = 0.0291          Pr(T > t) = 0.9855

.
. /* Clearly the SDs do differ somewhat, and allowing for this inequality using
> the unequal option increases the p-value a bit.
>
> We can also use a paired t-test to compare the LDL values measured at baseline
> to the year 1 value: */
.
. ttest ldl0 == ldl1

Paired t test
---------+--------------------------------------------------------------------
Variable |     Obs        Mean    Std. Err.   Std. Dev.   [95% Conf. Interval]
---------+--------------------------------------------------------------------
    ldl0 |   2,597     144.758    .7421193    37.81897    143.3028    146.2132
    ldl1 |   2,597    132.3248    .7661903    39.04565    130.8224    133.8272
---------+--------------------------------------------------------------------
    diff |   2,597    12.43327     .681349    34.72207    11.09723    13.76931
---------+--------------------------------------------------------------------
     mean(diff) = mean(ldl0 - ldl1)                               t =  18.2480
 Ho: mean(diff) = 0                              degrees of freedom =     2596

 Ha: mean(diff) < 0           Ha: mean(diff) != 0           Ha: mean(diff) > 0
 Pr(T < t) = 1.0000         Pr(|T| > |t|) = 0.0000          Pr(T > t) = 0.0000

.
. /* Note that the sample size goes down.  Why?
>
> We could also use a non-parametric Wilcoxon ranksum test to compare baseline
> LDL levels by smoking status; this method essentially operates on the ranks
> of LDL, rather than  directly on the measured LDL values.  For within-person
> comparisons of baseline and year 1 LDL, we can use the signrank test. */
.
. ranksum ldl0, by(smoking)

Two-sample Wilcoxon rank-sum (Mann-Whitney) test

     smoking |      obs    rank sum    expected
-------------+---------------------------------
          no |     2397     3265492   3299470.5
         yes |      355      522636    488657.5
-------------+---------------------------------
    combined |     2752     3788128     3788128

unadjusted variance   1.952e+08
adjustment for ties   -966.3953
                     ----------
adjusted variance     1.952e+08

Ho: ldl0(smoking==no) = ldl0(smoking==yes)
             z =  -2.432
    Prob > |z| =   0.0150

. signrank ldl0 = ldl1

Wilcoxon signed-rank test

        sign |      obs   sum ranks    expected
-------------+---------------------------------
    positive |     1714   2399964.5     1686741
    negative |      877    973517.5     1686741
        zero |        6          21          21
-------------+---------------------------------
         all |     2597     3373503     3373503

unadjusted variance   1.460e+09
adjustment for ties     -2366.5
adjustment for zeros     -22.75
                     ----------
adjusted variance     1.460e+09

Ho: ldl0 = ldl1
             z =  18.663
    Prob > |z| =   0.0000

.
. /* In both cases, the results are similar to the t-test results.  Another
> non-parametric approach is to compare the proportions that are greater than
> the overall median, using a chi-square test.  */
.
. median ldl0, by(smoking)

Median test

   Greater |
  than the |    current smoker
    median |        no        yes |     Total
-----------+----------------------+----------
        no |     1,222        157 |     1,379
       yes |     1,175        198 |     1,373
-----------+----------------------+----------
     Total |     2,397        355 |     2,752

          Pearson chi2(1) =   5.6437   Pr = 0.018

   Continuity corrected:
          Pearson chi2(1) =   5.3768   Pr = 0.020

.
. /* This command has options for how ties at the median are handled; check the
> help to learn about this. */
.
. help median

.
. /* For multi-category predictors, parametric and nonparametric alternatives
> include ANOVA and Kruskal-Wallis test respectively. */
.
. anova lpa0 raceth

                         Number of obs =      2,759    R-squared     =  0.0470
                         Root MSE      =    31.8468    Adj R-squared =  0.0463

                  Source | Partial SS         df         MS        F    Prob>F
              -----------+----------------------------------------------------
                   Model |  137855.61          2   68927.804     67.96  0.0000
                         |
                  raceth |  137855.61          2   68927.804     67.96  0.0000
                         |
                Residual |  2795180.5      2,756   1014.2164
              -----------+----------------------------------------------------
                   Total |  2933036.1      2,758   1063.4649

. kwallis lpa0, by(raceth)

Kruskal-Wallis equality-of-populations rank test

  +----------------------+-------+-----------+
  |               raceth |   Obs |  Rank Sum |
  +----------------------+-------+-----------+
  |                White | 2,447 |  3.28e+06 |
  |     African American |   218 | 417774.50 |
  | Latina, Asian, Other |    94 | 114079.00 |
  +----------------------+-------+-----------+

chi-squared =   109.553 with 2 d.f.
probability =     0.0001

chi-squared with ties =   109.653 with 2 d.f.
probability =     0.0001

.
. /* In this case, the two tests gives qualitatively identical assessments of
> the statistical significance of differences in Lp(a) by race/ethnicity.  Next
> week, we will look at the distribution of the these outcomes as well as sample
> size to help decide whether the parametric or non-parametric test is more
> appropriate. */
.
. log close
      name:  <unnamed>
       log:  /Users/eric/Work/MB/ATCR/200/Labs/lab 4.smcl
  log type:  smcl
 closed on:  20 Sep 2018, 11:29:52
-------------