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

. * do-file for lecture 1b of VHM 802/812, Winter 2014
. version 13

. set more off

. cd "f:\vhm\vhm802\data_stata"
f:\vhm\vhm802\data_stata

.
. use daisy2red, clear

. * multiple regression model with parity, twin, vag_disch
. regress milk120 parity twin vag_disch

      Source |       SS       df       MS              Number of obs =    1536
-------------+------------------------------           F(  3,  1532) =   90.76
       Model |   112904721     3  37634906.9           Prob > F      =  0.0000
    Residual |   635235472  1532  414644.564           R-squared     =  0.1509
-------------+------------------------------           Adj R-squared =  0.1493
       Total |   748140192  1535  487387.748           Root MSE      =  643.93

-------------+----------------------------------------------------------------
      milk120 |      Coef.    Std. Err.       t    P>|t|      [95% Con f. Interval]
-------------+----------------------------------------------------------------
      parity |     180.893    11.02613     16.41    0.000      159.2651     202.5209
        twin |   -197.0058    125.9608     -1.56    0.118     -444.0797     50.06806
   vag_disch |     199.717    74.45821      2.68    0.007      53.66624     345.7678
       _cons |    2713.174    34.68445     78.22    0.000       2645.14     2781.208
-------------+----------------------------------------------------------------

. * predictions, based on fitted values
. predict fit
(option xb assumed; fitted values)

. twoway (scatter milk120 parity, msymbol(circle_hollow) msize(vsmall)) (line fit parity)

. twoway (scatter milk120 parity, msymbol(circle_hollow) msize(vsmall)) (line fit parity, sort)

. * does not work because fitted values not on a single line (due to other predictors)
. twoway (scatter milk120 parity, msymbol(circle_hollow) msize(vsmall)) (lfit fit parity)

. twoway (scatter milk120 parity, msymbol(circle_hollow) msize(vsmall)) (lowess fit parity)

. * produces line/curve but is wrong (both commands)
. twoway (scatter milk120 parity, msymbol(circle_hollow) msize(vsmall)) (line fit parity, sort) if twin==0&vag_disch==0

. * ok but for a single group only
.
. * predictions, based on constructed fitted values
. set obs 1582 /* adding 8 rows to the worksheet */
obs was 1574, now 1582

. browse fit parity twin vag_disch in 1570/1582

. replace parity=1 in 1575/1578
(4 real changes made)

. replace parity=7 in 1579/1582
(4 real changes made)

. replace twin=0 in 1575/1582
(8 real changes made)

. replace twin=1 if mod(_n,2)==0 in 1575/1582
(4 real changes made)

. replace vag_disch=0 in 1575/1582
(8 real changes made)

. replace vag_disch=1 if _n==1577|_n==1578|_n==1581|_n==1582
(4 real changes made)

. drop fit

. predict fit /* fit now includes extra rows */
(option xb assumed; fitted values)

. twoway (line fit parity if twin==0&vag_disch==0 in 1575/1582) ///
>   (line fit parity if twin==0&vag_disch==1 in 1575/1582) ///
>   (line fit parity if twin==1&vag_disch==0 in 1575/1582) ///
>   (line fit parity if twin==1&vag_disch==1 in 1575/1582), legend(off)

.
. * model reduction
. regress milk120 parity /* gives SSE(R) and DFE(R) */

      Source |       SS       df       MS              Number of obs =    1536
-------------+------------------------------           F(  1,  1534) =  262.27
       Model |   109234227     1   109234227           Prob > F      =  0.0000
    Residual |   638905966  1534  416496.718           R-squared     =  0.1460
-------------+------------------------------           Adj R-squared =  0.1455
       Total |   748140192  1535  487387.748           Root MSE      =  645.37

-------------+----------------------------------------------------------------
      milk120 |      Coef.    Std. Err.       t    P>|t|      [95% Con f. Interval]
-------------+----------------------------------------------------------------
      parity |     178.347    11.01266     16.19    0.000      156.7455     199.9484
       _cons |     2727.08    34.33991     79.41    0.000      2659.722     2794.438
-------------+----------------------------------------------------------------

. regress milk120 parity twin vag_disch /* gives SSE(F) and DFE(F) */

      Source |       SS       df       MS              Number of obs =    1536
-------------+------------------------------           F(  3,  1532) =   90.76
       Model |   112904721     3  37634906.9           Prob > F      =  0.0000
    Residual |   635235472  1532  414644.564           R-squared     =  0.1509
-------------+------------------------------           Adj R-squared =  0.1493
       Total |   748140192  1535  487387.748           Root MSE      =  643.93

-------------+----------------------------------------------------------------
      milk120 |      Coef.    Std. Err.       t    P>|t|      [95% Con f. Interval]
-------------+----------------------------------------------------------------
      parity |     180.893    11.02613     16.41    0.000      159.2651     202.5209
        twin |   -197.0058    125.9608     -1.56    0.118     -444.0797     50.06806
   vag_disch |     199.717    74.45821      2.68    0.007      53.66624     345.7678
       _cons |    2713.174    34.68445     78.22    0.000       2645.14     2781.208
-------------+----------------------------------------------------------------

. scalar F=(638905966-635235472)/(1534-1532)/414644.564

. display "F=" F "  P=" Ftail(2,1532,F)
F=4.4260727  P=.01211471

. test twin vag_disch /* F-test for dropping twin and vag_disch */


 ( 1)  twin = 0


 ( 2)  vag_disch = 0


       F(  2,  1532) =    4.43
            Prob > F =    0.0121

.
. * model reduction, expanded example
. generate milk120k=milk120/1000
(46 missing values generated)

. regress milk120k parity twin vag_disch dyst rp

      Source |       SS       df       MS              Number of obs =    1536
-------------+------------------------------           F(  5,  1530) =   54.40
       Model |  112.932559     5  22.5865117           Prob > F      =  0.0000
    Residual |  635.207629  1530  .415168385           R-squared     =  0.1510
-------------+------------------------------           Adj R-squared =  0.1482
       Total |  748.140188  1535  .487387745           Root MSE      =  .64434

-------------+----------------------------------------------------------------
     milk120k |      Coef.    Std. Err.       t    P>|t|      [95% Con f. Interval]
-------------+----------------------------------------------------------------
      parity |     .180591    .0111486     16.20    0.000      .1587228     .2024591
        twin |   -.1938315    .1268007     -1.53    0.127     -.4425531       .05489
   vag_disch |    .2039539    .0763997      2.67    0.008      .0540946     .3538132
        dyst |   -.0116861    .0705298     -0.17    0.868     -.1500315     .1266592
          rp |   -.0110612    .0575734     -0.19    0.848     -.1239924     .1018699
       _cons |    2.715483    .0359492     75.54    0.000      2.644968     2.785998
-------------+----------------------------------------------------------------

. test dyst rp


 ( 1)  dyst = 0


 ( 2)  rp = 0


       F(  2,  1530) =    0.03
            Prob > F =    0.9670

. test twin dyst rp vag_disch


 ( 1)  twin = 0


 ( 2)  dyst = 0


 ( 3)  rp = 0


 ( 4)  vag_disch = 0


       F(  4,  1530) =    2.23
            Prob > F =    0.0640

. regress milk120k parity twin vag_disch

      Source |       SS       df       MS              Number of obs =    1536
-------------+------------------------------           F(  3,  1532) =   90.76
       Model |   112.90472     3  37.6349065           Prob > F      =  0.0000
    Residual |  635.235468  1532  .414644562           R-squared     =  0.1509
-------------+------------------------------           Adj R-squared =  0.1493
       Total |  748.140188  1535  .487387745           Root MSE      =  .64393

-------------+----------------------------------------------------------------
     milk120k |      Coef.    Std. Err.       t    P>|t|      [95% Con f. Interval]
-------------+----------------------------------------------------------------
      parity |     .180893    .0110261     16.41    0.000      .1592651     .2025209
        twin |   -.1970058    .1259608     -1.56    0.118     -.4440797      .050068
   vag_disch |     .199717    .0744582      2.68    0.007      .0536662     .3457678
       _cons |    2.713174    .0346844     78.22    0.000       2.64514     2.781208
-------------+----------------------------------------------------------------

. regress milk120k parity

      Source |       SS       df       MS              Number of obs =    1536
-------------+------------------------------           F(  1,  1534) =  262.27
       Model |  109.234225     1  109.234225           Prob > F      =  0.0000
    Residual |  638.905962  1534  .416496716           R-squared     =  0.1460
-------------+------------------------------           Adj R-squared =  0.1455
       Total |  748.140188  1535  .487387745           Root MSE      =  .64537

-------------+----------------------------------------------------------------
     milk120k |      Coef.    Std. Err.       t    P>|t|      [95% Con f. Interval]
-------------+----------------------------------------------------------------
      parity |     .178347    .0110127     16.19    0.000      .1567455     .1999484
       _cons |     2.72708    .0343399     79.41    0.000      2.659722     2.794438
-------------+----------------------------------------------------------------

. regress milk120k /* not really needed */

      Source |       SS       df       MS              Number of obs =    1536
-------------+------------------------------           F(  0,  1535) =    0.00
       Model |           0     0           .           Prob > F      =       .
    Residual |  748.140188  1535  .487387745           R-squared     =  0.0000
-------------+------------------------------           Adj R-squared =  0.0000
       Total |  748.140188  1535  .487387745           Root MSE      =  .69813

-------------+----------------------------------------------------------------
     milk120k |      Coef.    Std. Err.       t    P>|t|      [95% Con f. Interval]
-------------+----------------------------------------------------------------
       _cons |    3.215096    .0178132    180.49    0.000      3.180155     3.250036
-------------+----------------------------------------------------------------

.
. * polynomial regression
. generate parity2=parity^2

. regress milk120 parity parity2

      Source |       SS       df       MS              Number of obs =    1536
-------------+------------------------------           F(  2,  1533) =  212.78
       Model |   162559767     2  81279883.3           Prob > F      =  0.0000
    Residual |   585580426  1533  381983.318           R-squared     =  0.2173
-------------+------------------------------           Adj R-squared =  0.2163
       Total |   748140192  1535  487387.748           Root MSE      =  618.05

-------------+----------------------------------------------------------------
      milk120 |      Coef.    Std. Err.       t    P>|t|      [95% Con f. Interval]
-------------+----------------------------------------------------------------
      parity |    697.5025    45.18718     15.44    0.000      608.8673     786.1378
     parity2 |   -82.55675    6.987264    -11.82    0.000     -96.26236    -68.85115
       _cons |    2109.224    61.77411     34.14    0.000      1988.054     2230.395
-------------+----------------------------------------------------------------

. regress milk120 c.parity##c.parity

      Source |       SS       df       MS              Number of obs =    1536
-------------+------------------------------           F(  2,  1533) =  212.78
       Model |   162559767     2  81279883.3           Prob > F      =  0.0000
    Residual |   585580426  1533  381983.318           R-squared     =  0.2173
-------------+------------------------------           Adj R-squared =  0.2163
       Total |   748140192  1535  487387.748           Root MSE      =  618.05

------------------+----------------------------------------------------------------
           milk120 |      Coef.    Std. Err.       t    P>|t|      [95% Con f. Interval]
------------------+----------------------------------------------------------------
           parity |    697.5025    45.18718     15.44    0.000      608.8673     786.1378
                  |
c.parity#c.parity |   -82.55675    6.987264    -11.82    0.000     -96.26236    -68.85115
                  |
            _cons |    2109.224    61.77411     34.14    0.000      1988.054     2230.395
------------------+----------------------------------------------------------------

. * same model using factor notation, without generating parity2
.
. predict fit2
(option xb assumed; fitted values)

. twoway (scatter milk120 parity, msymbol(circle_hollow) msize(vsmall)) (line fit2 parity, sort)

. predict stdres, rstandard
(46 missing values generated)

. lowess stdres parity, xtitle("parity")

. generate parity3=parity^3

. regress milk120 parity parity2 parity3

      Source |       SS       df       MS              Number of obs =    1536
-------------+------------------------------           F(  3,  1532) =  163.18
       Model |   181175181     3  60391727.1           Prob > F      =  0.0000
    Residual |   566965011  1532    370081.6           R-squared     =  0.2422
-------------+------------------------------           Adj R-squared =  0.2407
       Total |   748140192  1535  487387.748           Root MSE      =  608.34

-------------+----------------------------------------------------------------
      milk120 |      Coef.    Std. Err.       t    P>|t|      [95% Con f. Interval]
-------------+----------------------------------------------------------------
      parity |    1635.566     139.543     11.72    0.000       1361.85     1909.281
     parity2 |   -409.1231    46.55594     -8.79    0.000     -500.4432     -317.803
     parity3 |    32.48599    4.580454      7.09    0.000      23.50137     41.47061
       _cons |    1396.606    117.4432     11.89    0.000      1166.239     1626.972
-------------+----------------------------------------------------------------

. predict fit3
(option xb assumed; fitted values)

. twoway (scatter milk120 parity, msymbol(circle_hollow) msize(vsmall)) (line fit3 parity, sort)

.
. * 1-way ANOVA
. oneway milk120k parity, tab

  Lactation |         Summary of milk120k
     number |        Mean   Std. Dev.       Freq.
------------+------------------------------------
          1 |   2.6396452   .48640371         403
          2 |   3.3478585   .62647033         369
          3 |   3.4294886   .65510671         308
          4 |   3.4881589   .65641493         219
          5 |   3.4272543   .69278456         164
          6 |   3.5178058   .54256333          69
          7 |      3.5656   .10381738           4
------------+------------------------------------
      Total |   3.2150956   .69813161        1536

                        Analysis of Variance
    Source              SS         df      MS            F     Prob > F
------------------------------------------------------------------------
Between groups      184.637128      6   30.7728546     83.50     0.0000
 Within groups       563.50306   1529   .368543532
------------------------------------------------------------------------
    Total           748.140188   1535   .487387745

Bartlett's test for equal variances:  chi2(6) =  55.9053  Prob>chi2 = 0.000

. regress milk120k parity

      Source |       SS       df       MS              Number of obs =    1536
-------------+------------------------------           F(  1,  1534) =  262.27
       Model |  109.234225     1  109.234225           Prob > F      =  0.0000
    Residual |  638.905962  1534  .416496716           R-squared     =  0.1460
-------------+------------------------------           Adj R-squared =  0.1455
       Total |  748.140188  1535  .487387745           Root MSE      =  .64537

-------------+----------------------------------------------------------------
     milk120k |      Coef.    Std. Err.       t    P>|t|      [95% Con f. Interval]
-------------+----------------------------------------------------------------
      parity |     .178347    .0110127     16.19    0.000      .1567455     .1999484
       _cons |     2.72708    .0343399     79.41    0.000      2.659722     2.794438
-------------+----------------------------------------------------------------

. scalar F_lackfit=(638.905962-563.50306)/(1534-1529)/0.368543532

. display "F=" F_lackfit "  P="Ftail(5,1529,F_lackfit)
F=40.919401  P=1.320e-39

.
. * simple example of collinearity in MER data bw5k
. use bw5k.dta, clear

. describe bwt cig*

              storage   display    value
variable name   type    format     label      variable label
-------------


bwt               double    %10.0g                      Birth weight (gms.)


cig_1             byte      %8.0g                       Cigarettes 1st trimester


cig_2             byte      %8.0g                       Cigarettes 2nd trimester


cig_3             byte      %8.0g                       Cigarettes 3rd trimester


. regress bwt cig_3

      Source |       SS       df       MS              Number of obs =    5000
-------------+------------------------------           F(  1,  4998) =   22.55
       Model |  7189497.98     1  7189497.98           Prob > F      =  0.0000
    Residual |  1.5935e+09  4998  318835.965           R-squared     =  0.0045
-------------+------------------------------           Adj R-squared =  0.0043
       Total |  1.6007e+09  4999  320210.372           Root MSE      =  564.66

-------------+----------------------------------------------------------------
          bwt |      Coef.    Std. Err.       t    P>|t|      [95% Con f. Interval]
-------------+----------------------------------------------------------------
       cig_3 |   -12.48752    2.629726     -4.75    0.000     -17.64293    -7.332101
       _cons |    3303.448    8.177729    403.96    0.000      3287.416      3319.48
-------------+----------------------------------------------------------------

. regress bwt cig_3 cig_1

      Source |       SS       df       MS              Number of obs =    5000
-------------+------------------------------           F(  2,  4997) =   12.43
       Model |  7922681.75     2  3961340.87           Prob > F      =  0.0000
    Residual |  1.5928e+09  4997  318753.046           R-squared     =  0.0049
-------------+------------------------------           Adj R-squared =  0.0046
       Total |  1.6007e+09  4999  320210.372           Root MSE      =  564.58

-------------+----------------------------------------------------------------
          bwt |      Coef.    Std. Err.       t    P>|t|      [95% Con f. Interval]
-------------+----------------------------------------------------------------
       cig_3 |   -6.481119    4.753749     -1.36    0.173     -15.80055     2.838316
       cig_1 |   -5.400689    3.560983     -1.52    0.129     -12.38178     1.580401
       _cons |    3304.627    8.213572    402.34    0.000      3288.525     3320.729
-------------+----------------------------------------------------------------

. regress bwt cig_3 cig_2

      Source |       SS       df       MS              Number of obs =    5000
-------------+------------------------------           F(  2,  4997) =   12.61
       Model |  8037757.94     2  4018878.97           Prob > F      =  0.0000
    Residual |  1.5927e+09  4997  318730.017           R-squared     =  0.0050
-------------+------------------------------           Adj R-squared =  0.0046
       Total |  1.6007e+09  4999  320210.372           Root MSE      =  564.56

-------------+----------------------------------------------------------------
          bwt |      Coef.    Std. Err.       t    P>|t|      [95% Con f. Interval]
-------------+----------------------------------------------------------------
       cig_3 |   -.0859205    8.043799     -0.01    0.991      -15.8553     15.68346
       cig_2 |   -12.11372    7.425483     -1.63    0.103     -26.67093      2.44348
       _cons |    3304.134    8.187191    403.57    0.000      3288.084     3320.185
-------------+----------------------------------------------------------------

. regress bwt cig_1 cig_2 cig_3

      Source |       SS       df       MS              Number of obs =    5000
-------------+------------------------------           F(  3,  4996) =    8.58
       Model |  8207559.77     3  2735853.26           Prob > F      =  0.0000
    Residual |  1.5925e+09  4996  318759.826           R-squared     =  0.0051
-------------+------------------------------           Adj R-squared =  0.0045
       Total |  1.6007e+09  4999  320210.372           Root MSE      =  564.59

-------------+----------------------------------------------------------------
          bwt |      Coef.    Std. Err.       t    P>|t|      [95% Con f. Interval]
-------------+----------------------------------------------------------------
       cig_1 |   -3.133501    4.293292     -0.73    0.466     -11.55024     5.283235
       cig_2 |   -8.463666     8.95284     -0.95    0.345     -26.01516     9.087831
       cig_3 |    -.337783    8.051574     -0.04    0.967      -16.1224     15.44684
       _cons |    3304.612    8.213676    402.33    0.000      3288.509     3320.714
-------------+----------------------------------------------------------------

. estat vif

    Variable |       VIF       1/VIF
-------------+----------------------
       cig_2 |     13.60    0.073506
       cig_3 |      9.38    0.106649
       cig_1 |      4.75    0.210477
-------------+----------------------
    Mean VIF |      9.24

. pwcorr cig*

             |    cig_1    cig_2    cig_3
-------------+---------------------------
       cig_1 |   1.0000
       cig_2 |   0.8883   1.0000
       cig_3 |   0.8331   0.9451   1.0000

. estat vce, corr

Correlation matrix of coefficients of regress model

e(V)   | cig_1    cig_2    cig_3    _cons
-------------+----------------------------------------
cig_1   |   1.0000
cig_2   |  -0.5586      1.0000
cig_3   |   0.0429     -0.8071      1.0000
_cons   |  -0.0797      0.0020     -0.0251      1.0000

.
. * more complex collinearity example in RC data Coleman
. use coleman.dta, clear

. summarize y x1-x5, sep(0)

    Variable |       Obs        Mean    Std. Dev.       Min        Max
-------------+--------------------------------------------------------
           y |        20     35.0825    5.817054       22.7       43.1
          x1 |        20      2.7315    .4540957       2.07       3.83
          x2 |        20      40.906    25.89855       9.99      86.27
          x3 |        20      3.1405    9.625382     -16.04      15.03
          x4 |        20      25.069    1.313805       21.6      28.01
          x5 |        20       6.255    .6543176       5.17       7.51

. pwcorr y x1-x5

             |        y       x1       x2       x3       x4       x5
-------------+------------------------------------------------------
           y |   1.0000
          x1 |   0.1923   1.0000
          x2 |   0.7534   0.1811   1.0000
          x3 |   0.9272   0.2296   0.8272   1.0000
          x4 |   0.3336   0.5027   0.0511   0.1833   1.0000
          x5 |   0.7330   0.1968   0.9271   0.8191   0.1238   1.0000

. regress y x1-x5

      Source |       SS       df       MS              Number of obs =      20
-------------+------------------------------           F(  5,    14) =   27.08
       Model |  582.686401     5   116.53728           Prob > F      =  0.0000
    Residual |  60.2378924    14   4.3027066           R-squared     =  0.9063
-------------+------------------------------           Adj R-squared =  0.8728
       Total |  642.924294    19  33.8381207           Root MSE      =  2.0743

-------------+----------------------------------------------------------------
            y |      Coef.    Std. Err.       t    P>|t|      [95% Con f. Interval]
-------------+----------------------------------------------------------------
          x1 |   -1.793333    1.233397     -1.45    0.168     -4.438706     .8520392
          x2 |    .0436015    .0532589      0.82    0.427     -.0706275     .1578304
          x3 |    .5557601    .0929564      5.98    0.000      .3563884     .7551318
          x4 |    1.110168    .4337681      2.56    0.023      .1798279     2.040508
          x5 |   -1.810919     2.02739     -0.89    0.387     -6.159239       2.5374
       _cons |    19.94857    13.62755      1.46    0.165     -9.279627     49.17676
-------------+----------------------------------------------------------------

. estat vif

    Variable |       VIF       1/VIF
-------------+----------------------
          x2 |      8.40    0.119029
          x5 |      7.77    0.128687
          x3 |      3.54    0.282874
          x4 |      1.43    0.697286
          x1 |      1.39    0.721918
-------------+----------------------
    Mean VIF |      4.51

. estat vce, corr

Correlation matrix of coefficients of regress model

e(V)   | x1    x2    x3    x4    x5    _cons
-------------+------------------------------------------------------------
x1   |   1.0000
x2   |  -0.0785      1.0000
x3   |  -0.0282     -0.3454      1.0000
x4   |  -0.4860      0.2392     -0.1713      1.0000
x5   |   0.0340     -0.7794     -0.2080     -0.1488      1.0000
_cons   |   0.1221      0.4013      0.3710     -0.5739     -0.6912      1.0000

. * note: code for reduced models left for the interested student
.
end of do-file

. exit, clear