. do "C:\Users\DEFAUL~1.SID\AppData\Local\Temp\STD00000000.tmp"

. * do-file for lecture 6 of VHM 802, Winter 2016
. version 14 /* works also with version 13 */

. set more off

. cd "h:\vhm\vhm802\data_csv"
h:\vhm\vhm802\data_csv

.
. import delimited ch03ex1.csv, clear
(2 vars, 29 obs)

. oneway liverwgt diet, tabulate

            |         Summary of liverwgt
       diet |        Mean   Std. Dev.       Freq.
------------+------------------------------------
          1 |   3.7457143   .28401046           7
          2 |        3.58   .18213024           8
          3 |   3.5983333   .09621159           6
          4 |      3.9225   .19710402           8
------------+------------------------------------
      Total |   3.7182759   .23998614          29

                        Analysis of Variance
    Source              SS         df      MS            F     Prob > F
------------------------------------------------------------------------
Between groups      .578208876      3   .192736292      4.66     0.0102
 Within groups       1.0344049     25   .041376196
------------------------------------------------------------------------
    Total           1.61261378     28   .057593349

Bartlett's test for equal variances:  chi2(3) =   5.1211  Prob>chi2 = 0.163

. anova liverwgt diet /* allows postestimation commands */

                         Number of obs =         29    R-squared     =  0.3586
                         Root MSE      =    .203411    Adj R-squared =  0.2816

                  Source | Partial SS         df         MS        F    Prob>F
              -----------+----------------------------------------------------
                   Model |  .57820888          3   .19273629      4.66  0.0102
                         |
                    diet |  .57820888          3   .19273629      4.66  0.0102
                         |
                Residual |  1.0344049         25    .0413762
              -----------+----------------------------------------------------
                   Total |  1.6126138         28   .05759335

. regress /* estimates corresponding to anova model */

      Source |       SS           df       MS      Number of obs   =        29
-------------+----------------------------------   F(3, 25)        =      4.66
       Model |  .578208876         3  .192736292   Prob > F        =    0.0102
    Residual |   1.0344049        25  .041376196   R-squared       =    0.3586
-------------+----------------------------------   Adj R-squared   =    0.2816
       Total |  1.61261378        28  .057593349   Root MSE        =    .20341

-------------+----------------------------------------------------------------
     liverwgt |      Coef.    Std. Err.       t    P>|t|      [95% Con f. Interval]
-------------+----------------------------------------------------------------
        diet |
          2  |   -.1657143    .1052754     -1.57    0.128      -.382533     .0511045
          3  |    -.147381    .1131677     -1.30    0.205     -.3804541     .0856922
          4  |    .1767857    .1052754      1.68    0.106     -.0400331     .3936044
             |
       _cons |    3.745714    .0768823     48.72    0.000      3.587372     3.904056
-------------+----------------------------------------------------------------

. regress liverwgt i.diet /* totally identical */

      Source |       SS           df       MS      Number of obs   =        29
-------------+----------------------------------   F(3, 25)        =      4.66
       Model |  .578208876         3  .192736292   Prob > F        =    0.0102
    Residual |   1.0344049        25  .041376196   R-squared       =    0.3586
-------------+----------------------------------   Adj R-squared   =    0.2816
       Total |  1.61261378        28  .057593349   Root MSE        =    .20341

-------------+----------------------------------------------------------------
     liverwgt |      Coef.    Std. Err.       t    P>|t|      [95% Con f. Interval]
-------------+----------------------------------------------------------------
        diet |
          2  |   -.1657143    .1052754     -1.57    0.128      -.382533     .0511045
          3  |    -.147381    .1131677     -1.30    0.205     -.3804541     .0856922
          4  |    .1767857    .1052754      1.68    0.106     -.0400331     .3936044
             |
       _cons |    3.745714    .0768823     48.72    0.000      3.587372     3.904056
-------------+----------------------------------------------------------------

. xi: boxcox liverwgt i.diet /* Box-Cox analysis; note: needs xi: */
i.diet            _Idiet_1-4          (naturally coded; _Idiet_1 omitted)
Fitting comparison model

Iteration 0:   log likelihood =  .74765532
Iteration 1:   log likelihood =  1.6767683
Iteration 2:   log likelihood =  1.6782923
Iteration 3:   log likelihood =  1.6782967
Iteration 4:   log likelihood =  1.6782967

Fitting full model

Iteration 0:   log likelihood =  7.1860909
Iteration 1:   log likelihood =  8.0100122
Iteration 2:   log likelihood =  8.0113804
Iteration 3:   log likelihood =  8.0113804

                                                  Number of obs   =          29
                                                  LR chi2(3)      =      12.67
Log likelihood = 8.0113804                        Prob > chi2     =      0.005

-------------+----------------------------------------------------------------
    liverwgt |       Coef.   Std. Err.      z    P>|z|    [95% Conf. Interval]
-------------+----------------------------------------------------------------
      /theta |  -2.161941   2.514124    -0.86   0.390    -7.089534    2.765652
-------------+----------------------------------------------------------------

Estimates of scale-variant parameters
-------------+--------------
             |      Coef.
-------------+--------------
Notrans      |
    _Idiet_2 |  -.0025137
    _Idiet_3 |  -.0020231
    _Idiet_4 |    .002801
       _cons |   .4354793
-------------+--------------
      /sigma |   .0029048
-------------+--------------

---------------------------------------------------------
   Test         Restricted     LR statistic      P-value
    H0:       log likelihood       chi2       Prob > chi2
---------------------------------------------------------
theta = -1       7.9027094         0.22           0.641
theta =  0       7.6300907         0.76           0.383
theta =  1       7.1860909         1.65           0.199
---------------------------------------------------------

. * means and SE with CI
. anova liverwgt diet

                         Number of obs =         29    R-squared     =  0.3586
                         Root MSE      =    .203411    Adj R-squared =  0.2816

                  Source | Partial SS         df         MS        F    Prob>F
              -----------+----------------------------------------------------
                   Model |  .57820888          3   .19273629      4.66  0.0102
                         |
                    diet |  .57820888          3   .19273629      4.66  0.0102
                         |
                Residual |  1.0344049         25    .0413762
              -----------+----------------------------------------------------
                   Total |  1.6126138         28   .05759335

. lincom _cons+1.diet /* similar for diets 2-4 */


 ( 1)  1b.diet + _cons = 0


-------------+----------------------------------------------------------------
     liverwgt |      Coef.    Std. Err.       t    P>|t|      [95% Con f. Interval]
-------------+----------------------------------------------------------------
         (1) |    3.745714    .0768823     48.72    0.000      3.587372     3.904056
-------------+----------------------------------------------------------------

. margins diet

Adjusted predictions                            Number of obs     =         29


Expression  : Linear prediction, predict()


-------------+----------------------------------------------------------------
             |            Delta-method
             |     Margin    Std. Err.       t    P>|t|      [95% Con f. Interval]
-------------+----------------------------------------------------------------
        diet |
          1  |    3.745714    .0768823     48.72    0.000      3.587372     3.904056
          2  |        3.58    .0719168     49.78    0.000      3.431885     3.728115
          3  |    3.598333    .0830424     43.33    0.000      3.427304     3.769362
          4  |      3.9225    .0719168     54.54    0.000      3.774385     4.070615
-------------+----------------------------------------------------------------

. marginsplot /* interval plot */


Variables that uniquely identify margins: diet


. * pairwise comparisons
. oneway liverwgt diet, bonferroni /* bonferroni comparisons */

                        Analysis of Variance
    Source              SS         df      MS            F     Prob > F
------------------------------------------------------------------------
Between groups      .578208876      3   .192736292      4.66     0.0102
 Within groups       1.0344049     25   .041376196
------------------------------------------------------------------------
    Total           1.61261378     28   .057593349

Bartlett's test for equal variances:  chi2(3) =   5.1211  Prob>chi2 = 0.163

                        Comparison of liverwgt by diet
                                (Bonferroni)
Row Mean-|
Col Mean |          1          2          3
---------+---------------------------------
       2 |   -.165714
         |      0.768
         |
       3 |   -.147381    .018333
         |      1.000      1.000
         |
       4 |    .176786      .3425    .324167
         |      0.633      0.015      0.041

. * general method for anova and regression
. anova liverwgt diet

                         Number of obs =         29    R-squared     =  0.3586
                         Root MSE      =    .203411    Adj R-squared =  0.2816

                  Source | Partial SS         df         MS        F    Prob>F
              -----------+----------------------------------------------------
                   Model |  .57820888          3   .19273629      4.66  0.0102
                         |
                    diet |  .57820888          3   .19273629      4.66  0.0102
                         |
                Residual |  1.0344049         25    .0413762
              -----------+----------------------------------------------------
                   Total |  1.6126138         28   .05759335

. pwcompare diet, pv mcomp(noadjust) /* no adjustment - the default */

Pairwise comparisons of marginal linear predictions


Margins  : asbalanced


-------------+---------------------------------------
             |                            Una djusted
             |   Contrast    Std. Err.       t    P>|t|
-------------+---------------------------------------
        diet |
     2 vs 1  |   -.1657143    .1052754     -1.57    0.128
     3 vs 1  |    -.147381    .1131677     -1.30    0.205
     4 vs 1  |    .1767857    .1052754      1.68    0.106
     3 vs 2  |    .0183333    .1098547      0.17    0.869
     4 vs 2  |    .3424999    .1017057      3.37    0.002
     4 vs 3  |    .3241666    .1098547      2.95    0.007
-------------+---------------------------------------

. pwcompare diet, pv mcomp(bon)      /* Bonferroni methodt */

Pairwise comparisons of marginal linear predictions


Margins  : asbalanced


-------------+-------------
             |    Number of
             |  Comparisons
-------------+-------------
        diet |             6
-------------+-------------

-------------+---------------------------------------
             |                            Bon ferroni
             |   Contrast    Std. Err.       t    P>|t|
-------------+---------------------------------------
        diet |
     2 vs 1  |   -.1657143    .1052754     -1.57    0.768
     3 vs 1  |    -.147381    .1131677     -1.30    1.000
     4 vs 1  |    .1767857    .1052754      1.68    0.633
     3 vs 2  |    .0183333    .1098547      0.17    1.000
     4 vs 2  |    .3424999    .1017057      3.37    0.015
     4 vs 3  |    .3241666    .1098547      2.95    0.041
-------------+---------------------------------------

. pwcompare diet, pv mcomp(tukey)    /* Tukey method */

Pairwise comparisons of marginal linear predictions


Margins  : asbalanced


-------------+-------------
             |    Number of
             |  Comparisons
-------------+-------------
        diet |             6
-------------+-------------

-------------+---------------------------------------
             |                              T ukey
             |   Contrast    Std. Err.       t    P>|t|
-------------+---------------------------------------
        diet |
     2 vs 1  |   -.1657143    .1052754     -1.57    0.411
     3 vs 1  |    -.147381    .1131677     -1.30    0.570
     4 vs 1  |    .1767857    .1052754      1.68    0.355
     3 vs 2  |    .0183333    .1098547      0.17    0.998
     4 vs 2  |    .3424999    .1017057      3.37    0.012
     4 vs 3  |    .3241666    .1098547      2.95    0.032
-------------+---------------------------------------

. * general approach to testing
. anova liverwgt diet

                         Number of obs =         29    R-squared     =  0.3586
                         Root MSE      =    .203411    Adj R-squared =  0.2816

                  Source | Partial SS         df         MS        F    Prob>F
              -----------+----------------------------------------------------
                   Model |  .57820888          3   .19273629      4.66  0.0102
                         |
                    diet |  .57820888          3   .19273629      4.66  0.0102
                         |
                Residual |  1.0344049         25    .0413762
              -----------+----------------------------------------------------
                   Total |  1.6126138         28   .05759335

. test, showorder

 Order of columns in the design matrix
1:       (diet==1)
2:       (diet==2)
3:       (diet==3)
4:       (diet==4)
5:       _cons

. matrix input mycon=(1,-1,0,0,0\1,0,-1,0,0\1,0,0,-1,0\0,1,-1,0,0\0,1,0,-1,0\0,0,1,-1,0)

. test, test(mycon) mtest


 ( 1)  1b.diet - 2.diet = 0


 ( 2)  1b.diet - 3.diet = 0


 ( 3)  1b.diet - 4.diet = 0


 ( 4)  2.diet - 3.diet = 0


 ( 5)  2.diet - 4.diet = 0


 ( 6)  3.diet - 4.diet = 0


       Constraint 3 dropped
       Constraint 4 dropped
       Constraint 6 dropped

-------+-------------------------------
       |F(df,25)       df       p
-------+-------------------------------
  (1)  |        2.48      1     0.1280 #
  (2)  |        1.70      1     0.2047 #
  (3)  |        2.82      1     0.1056 #
  (4)  |        0.03      1     0.8688 #
  (5)  |       11.34      1     0.0025 #
  (6)  |        8.71      1     0.0068 #
-------+-------------------------------
  all  |        4.66      3     0.0102
-------+-------------------------------
# unadjusted p-values

. test, test(mycon) mtest(bon)  /* Bonferroni method*/


 ( 1)  1b.diet - 2.diet = 0


 ( 2)  1b.diet - 3.diet = 0


 ( 3)  1b.diet - 4.diet = 0


 ( 4)  2.diet - 3.diet = 0


 ( 5)  2.diet - 4.diet = 0


 ( 6)  3.diet - 4.diet = 0


       Constraint 3 dropped
       Constraint 4 dropped
       Constraint 6 dropped

-------+-------------------------------
       |F(df,25)       df       p
-------+-------------------------------
  (1)  |        2.48      1     0.7682 #
  (2)  |        1.70      1     1.0000 #
  (3)  |        2.82      1     0.6333 #
  (4)  |        0.03      1     1.0000 #
  (5)  |       11.34      1     0.0147 #
  (6)  |        8.71      1     0.0408 #
-------+-------------------------------
  all  |        4.66      3     0.0102
-------+-------------------------------
# Bonferroni-adjusted p-values

. test, test(mycon) mtest(holm) /* Holm method*/


 ( 1)  1b.diet - 2.diet = 0


 ( 2)  1b.diet - 3.diet = 0


 ( 3)  1b.diet - 4.diet = 0


 ( 4)  2.diet - 3.diet = 0


 ( 5)  2.diet - 4.diet = 0


 ( 6)  3.diet - 4.diet = 0


       Constraint 3 dropped
       Constraint 4 dropped
       Constraint 6 dropped

-------+-------------------------------
       |F(df,25)       df       p
-------+-------------------------------
  (1)  |        2.48      1     0.3841 #
  (2)  |        1.70      1     0.4094 #
  (3)  |        2.82      1     0.4222 #
  (4)  |        0.03      1     0.8688 #
  (5)  |       11.34      1     0.0147 #
  (6)  |        8.71      1     0.0340 #
-------+-------------------------------
  all  |        4.66      3     0.0102
-------+-------------------------------
# Holm-adjusted p-values

. * contrasts
. lincom (1.diet+2.diet+3.diet)/3-4.diet


 ( 1)  .3333333*1b.diet + .3333333*2.diet + .3333333*3.diet - 4.diet = 0


-------------+----------------------------------------------------------------
     liverwgt |      Coef.    Std. Err.       t    P>|t|      [95% Con f. Interval]
-------------+----------------------------------------------------------------
         (1) |   -.2811507     .084674     -3.32    0.003     -.4555401    -.1067614
-------------+----------------------------------------------------------------

. scalar tval=-.2811507/.084674 /* to get more decimals than in listing */

. di "SS: " tval^2*0.0413762 "  in %: " tval^2*0.0413762/.57820888*100
SS: .45617217  in %: 78.89401

. di "Scheffe test:  F = " tval^2/3 "  P = " Ftail(3,25,tval^2/3)
Scheffe test:  F = 3.6749965  P = .02546919

. lincom (2.diet+3.diet)/2-1.diet


 ( 1)  - 1b.diet + .5*2.diet + .5*3.diet = 0


-------------+----------------------------------------------------------------
     liverwgt |      Coef.    Std. Err.       t    P>|t|      [95% Con f. Interval]
-------------+----------------------------------------------------------------
         (1) |   -.1565476    .0944876     -1.66    0.110     -.3511484     .0380532
-------------+----------------------------------------------------------------

. lincom 2.diet-3.diet


 ( 1)  2.diet - 3.diet = 0


-------------+----------------------------------------------------------------
     liverwgt |      Coef.    Std. Err.       t    P>|t|      [95% Con f. Interval]
-------------+----------------------------------------------------------------
         (1) |   -.0183333    .1098547     -0.17    0.869     -.2445833     .2079167
-------------+----------------------------------------------------------------

.
. import delimited ch03ta1.csv, clear
(2 vars, 37 obs)

. regress logtime temp

      Source |       SS           df       MS      Number of obs   =        37
-------------+----------------------------------   F(1, 35)        =    325.41
       Model |  3.45926188         1  3.45926188   Prob > F        =    0.0000
    Residual |  .372062445        35  .010630356   R-squared       =    0.9029
-------------+----------------------------------   Adj R-squared   =    0.9001
       Total |  3.83132433        36  .106425676   Root MSE        =     .1031

-------------+----------------------------------------------------------------
      logtime |      Coef.    Std. Err.       t    P>|t|      [95% Con f. Interval]
-------------+----------------------------------------------------------------
        temp |   -.0118567    .0006573    -18.04    0.000     -.0131911    -.0105224
       _cons |    3.956007    .1391174     28.44    0.000      3.673584     4.238431
-------------+----------------------------------------------------------------

. anova logtime temp

                         Number of obs =         37    R-squared     =  0.9233
                         Root MSE      =    .095801    Adj R-squared =  0.9138

                  Source | Partial SS         df         MS        F    Prob>F
              -----------+----------------------------------------------------
                   Model |   3.537632          4   .88440801     96.36  0.0000
                         |
                    temp |   3.537632          4   .88440801     96.36  0.0000
                         |
                Residual |  .29369229         32   .00917788
              -----------+----------------------------------------------------
                   Total |  3.8313243         36   .10642568

. scalar F=(.3721-.2937)/(35-32)/.009178 /* F-value for lack of fit test */

. display F _newline Ftail(3,32,F) /* display F and P-value */
2.8473887
.05297014

. * polynomial contrasts
. lincom -2*175.temp-1*194.temp+0*213.temp+1*231.temp+2*250.temp /* linear */


 ( 1)  - 2*175b.temp - 194.temp + 231.temp + 2*250.temp = 0


-------------+----------------------------------------------------------------
      logtime |      Coef.    Std. Err.       t    P>|t|      [95% Con f. Interval]
-------------+----------------------------------------------------------------
         (1) |   -2.186131    .1147427    -19.05    0.000     -2.419854    -1.952408
-------------+----------------------------------------------------------------

. lincom 2*175.temp-1*194.temp-2*213.temp-1*231.temp+2*250.temp /* quadratic */


 ( 1)  2*175b.temp - 194.temp - 2*213.temp - 231.temp + 2*250.temp = 0


-------------+----------------------------------------------------------------
      logtime |      Coef.    Std. Err.       t    P>|t|      [95% Con f. Interval]
-------------+----------------------------------------------------------------
         (1) |    .4002977    .1332473      3.00    0.005      .1288819     .6717135
-------------+----------------------------------------------------------------

. lincom -1*175.temp+2*194.temp+0*213.temp-2*231.temp+1*250.temp /* cubic */


 ( 1)  - 175b.temp + 2*194.temp - 2*231.temp + 250.temp = 0


-------------+----------------------------------------------------------------
      logtime |      Coef.    Std. Err.       t    P>|t|      [95% Con f. Interval]
-------------+----------------------------------------------------------------
         (1) |   -.0069048    .1118496     -0.06    0.951     -.2347349     .2209254
-------------+----------------------------------------------------------------

. lincom 1*175.temp-4*194.temp+6*213.temp-4*231.temp+1*250.temp /* quartic */


 ( 1)  175b.temp - 4*194.temp + 6*213.temp - 4*231.temp + 250.temp = 0


-------------+----------------------------------------------------------------
      logtime |      Coef.    Std. Err.       t    P>|t|      [95% Con f. Interval]
-------------+----------------------------------------------------------------
         (1) |   -.0379762    .2886367     -0.13    0.896     -.6259099     .5499576
-------------+----------------------------------------------------------------

.
end of do-file

. exit, clear