. do "C:\Users\hstryhn\AppData\Local\Temp\STD26a8_000000.tmp"

. * do-file for lecture 2a of VHM 812/802, Winter 2025
. version 18 /* works mostly with versions 14-17, as indicated below */

. set more off

. set scheme stcolor_alt

. cd "r:\"
r:\

.
. use daisy2red, clear

. * recoding a continuous predictor (VER 14.4)
. regress milk120 parity

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

-------------+----------------------------------------------------------------
      milk120 | Coefficient   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
-------------+----------------------------------------------------------------

. generate parity1=parity-1

. regress milk120 parity1 // intercept now for parity=1

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

-------------+----------------------------------------------------------------
      milk120 | Coefficient   Std. err.       t    P>|t|      [95% con f. interval]
-------------+----------------------------------------------------------------
     parity1 |     178.347    11.01266     16.19    0.000      156.7455     199.9484
       _cons |    2905.427    25.23474    115.14    0.000      2855.928     2954.925
-------------+----------------------------------------------------------------

.
. * switch now to wpc outcome
. summarize wpc, d

           Interval from wait period to conception
-------------------------------------------------------------
      Percentiles      Smallest
 1%           11              1
 5%           15              1
10%           20              3       Obs               1,574
25%           30              4       Sum of wgt.       1,574

50%           53                      Mean           68.79924
                        Largest       Std. dev.      51.59928
75%           94            254
90%          141            255       Variance       2662.486
95%          178            284       Skewness       1.356612
99%          232            298       Kurtosis       4.551503

. * scaling of continuous predictor
. regress wpc herd_size

      Source |       SS           df       MS      Number of obs   =     1,574
-------------+----------------------------------   F(1, 1572)      =     70.66
       Model |  180151.608         1  180151.608   Prob > F        =    0.0000
    Residual |  4007938.95     1,572  2549.57949   R-squared       =    0.0430
-------------+----------------------------------   Adj R-squared   =    0.0424
       Total |  4188090.56     1,573  2662.48605   Root MSE        =    50.493

-------------+----------------------------------------------------------------
          wpc | Coefficient   Std. err.       t    P>|t|      [95% con f. interval]
-------------+----------------------------------------------------------------
   herd_size |    .1725618    .0205286      8.41    0.000      .1322954     .2128281
       _cons |    25.48492     5.30769      4.80    0.000      15.07403     35.89582
-------------+----------------------------------------------------------------

. generate hs100=herd_size/100

. regress wpc hs100 // slope corresponds to change in herd size of 100

      Source |       SS           df       MS      Number of obs   =     1,574
-------------+----------------------------------   F(1, 1572)      =     70.66
       Model |  180151.593         1  180151.593   Prob > F        =    0.0000
    Residual |  4007938.97     1,572   2549.5795   R-squared       =    0.0430
-------------+----------------------------------   Adj R-squared   =    0.0424
       Total |  4188090.56     1,573  2662.48605   Root MSE        =    50.493

-------------+----------------------------------------------------------------
          wpc | Coefficient   Std. err.       t    P>|t|      [95% con f. interval]
-------------+----------------------------------------------------------------
       hs100 |    17.25617    2.052862      8.41    0.000      13.22954     21.28281
       _cons |    25.48492     5.30769      4.80    0.000      15.07403     35.89582
-------------+----------------------------------------------------------------

. gen hs100_sq=hs100^2

. regress wpc hs100 hs100_sq // VER 14.8

      Source |       SS           df       MS      Number of obs   =     1,574
-------------+----------------------------------   F(2, 1571)      =     40.45
       Model |  205120.743         2  102560.371   Prob > F        =    0.0000
    Residual |  3982969.82     1,571   2535.3086   R-squared       =    0.0490
-------------+----------------------------------   Adj R-squared   =    0.0478
       Total |  4188090.56     1,573  2662.48605   Root MSE        =    50.352

-------------+----------------------------------------------------------------
          wpc | Coefficient   Std. err.       t    P>|t|      [95% con f. interval]
-------------+----------------------------------------------------------------
       hs100 |   -29.51616    15.04392     -1.96    0.050     -59.02444    -.0078794
    hs100_sq |    9.744227    3.104995      3.14    0.002      3.653856      15.8346
       _cons |    77.74837    17.47457      4.45    0.000      43.47243     112.0243
-------------+----------------------------------------------------------------

. regress wpc c.hs100##c.hs100 // same with Stata factor notation, see 2aL-14

      Source |       SS           df       MS      Number of obs   =     1,574
-------------+----------------------------------   F(2, 1571)      =     40.45
       Model |  205120.743         2  102560.372   Prob > F        =    0.0000
    Residual |  3982969.82     1,571   2535.3086   R-squared       =    0.0490
-------------+----------------------------------   Adj R-squared   =    0.0478
       Total |  4188090.56     1,573  2662.48605   Root MSE        =    50.352

----------------+----------------------------------------------------------------
             wpc | Coefficient   Std. err.       t    P>|t|      [95% con f. interval]
----------------+----------------------------------------------------------------
          hs100 |   -29.51616    15.04392     -1.96    0.050     -59.02444    -.0078797
                |
c.hs100#c.hs100 |    9.744227    3.104995      3.14    0.002      3.653856      15.8346
                |
          _cons |    77.74837    17.47457      4.45    0.000      43.47243     112.0243
----------------+----------------------------------------------------------------

. * no exploration of centering hs100 because not really needed
.
. * parametrizing categorical predictors, back to milk120 again...
. regress milk120 i.parity // VER 14.7

      Source |       SS           df       MS      Number of obs   =     1,536
-------------+----------------------------------   F(6, 1529)      =     83.50
       Model |   184637128         6  30772854.7   Prob > F        =    0.0000
    Residual |   563503064     1,529  368543.534   R-squared       =    0.2468
-------------+----------------------------------   Adj R-squared   =    0.2438
       Total |   748140192     1,535  487387.748   Root MSE        =    607.08

-------------+----------------------------------------------------------------
      milk120 | Coefficient   Std. err.       t    P>|t|      [95% con f. interval]
-------------+----------------------------------------------------------------
      parity |
          2  |    708.2134    43.74086     16.19    0.000      622.4149     794.0118
          3  |    789.8435    45.94638     17.19    0.000      699.7189     879.9681
          4  |    848.5137    50.96418     16.65    0.000      748.5467     948.4808
          5  |    787.6091    56.22915     14.01    0.000      677.3147     897.9035
          6  |    878.1606    79.09305     11.10    0.000      723.0183     1033.303
          7  |    925.9548    305.0416      3.04    0.002      327.6106     1524.299
             |
       _cons |    2639.645     30.2407     87.29    0.000      2580.328     2698.963
-------------+----------------------------------------------------------------

. tabstat milk120, by(parity)

Summary for variables: milk120
Group variable: parity (Lactation number)

parity   |      Mean
---------+----------
1   |  2639.645
2   |  3347.859
3   |  3429.489
4   |  3488.159
5   |  3427.254
6   |  3517.806
7   |    3565.6
---------+----------
Total   |  3215.096
---------+----------

. testparm i.parity // overall test for parity


 ( 1)  2.parity = 0


 ( 2)  3.parity = 0


 ( 3)  4.parity = 0


 ( 4)  5.parity = 0


 ( 5)  6.parity = 0


 ( 6)  7.parity = 0


       F(  6,  1529) =   83.50
            Prob > F =    0.0000

. regress milk120 ib7.parity // parametrized with last category (7) as reference

      Source |       SS           df       MS      Number of obs   =     1,536
-------------+----------------------------------   F(6, 1529)      =     83.50
       Model |   184637128         6  30772854.7   Prob > F        =    0.0000
    Residual |   563503064     1,529  368543.534   R-squared       =    0.2468
-------------+----------------------------------   Adj R-squared   =    0.2438
       Total |   748140192     1,535  487387.748   Root MSE        =    607.08

-------------+----------------------------------------------------------------
      milk120 | Coefficient   Std. err.       t    P>|t|      [95% con f. interval]
-------------+----------------------------------------------------------------
      parity |
          1  |   -925.9548    305.0416     -3.04    0.002     -1524.299    -327.6106
          2  |   -217.7414    305.1797     -0.71    0.476     -816.3565     380.8736
          3  |   -136.1113    305.5036     -0.45    0.656     -735.3618     463.1391
          4  |   -77.44108    306.2984     -0.25    0.800     -678.2506     523.3684
          5  |   -138.3457    307.2183     -0.45    0.653     -740.9596     464.2682
          6  |   -47.79418    312.2132     -0.15    0.878     -660.2056     564.6173
             |
       _cons |      3565.6    303.5389     11.75    0.000      2970.203     4160.997
-------------+----------------------------------------------------------------

. testparm i.parity // the same, of course


 ( 1)  1.parity = 0


 ( 2)  2.parity = 0


 ( 3)  3.parity = 0


 ( 4)  4.parity = 0


 ( 5)  5.parity = 0


 ( 6)  6.parity = 0


       F(  6,  1529) =   83.50
            Prob > F =    0.0000

.
. * exploring additivity and interaction
. * two categorical predictors (VER 14.9)
. regress wpc i.rp i.vag_disch // additive model

      Source |       SS           df       MS      Number of obs   =     1,574
-------------+----------------------------------   F(2, 1571)      =      4.65
       Model |  24663.1482         2  12331.5741   Prob > F        =    0.0097
    Residual |  4163427.41     1,571  2650.17658   R-squared       =    0.0059
-------------+----------------------------------   Adj R-squared   =    0.0046
       Total |  4188090.56     1,573  2662.48605   Root MSE        =     51.48

-------------+----------------------------------------------------------------
          wpc | Coefficient   Std. err.       t    P>|t|      [95% con f. interval]
-------------+----------------------------------------------------------------
          rp |
        yes  |     10.2397    4.540774      2.26    0.024      1.333087     19.14632
             |
   vag_disch |
        yes  |    9.066942      5.9819      1.52    0.130     -2.666406     20.80029
       _cons |    67.35756    1.381095     48.77    0.000      64.64857     70.06654
-------------+----------------------------------------------------------------

. regress wpc i.rp i.vag_disch i.rp#i.vag_disch // interaction added

      Source |       SS           df       MS      Number of obs   =     1,574
-------------+----------------------------------   F(3, 1570)      =      4.53
       Model |  35915.9774         3  11971.9925   Prob > F        =    0.0036
    Residual |  4152174.58     1,570  2644.69719   R-squared       =    0.0086
-------------+----------------------------------   Adj R-squared   =    0.0067
       Total |  4188090.56     1,573  2662.48605   Root MSE        =    51.427

-------------+----------------------------------------------------------------
          wpc | Coefficient   Std. err.       t    P>|t|      [95% con f. interval]
-------------+----------------------------------------------------------------
          rp |
        yes  |    6.339794    4.914322      1.29    0.197     -3.299531     15.97912
             |
   vag_disch |
        yes  |    .5429296    7.265382      0.07    0.940     -13.70794      14.7938
             |
rp#vag_disch |
    yes#yes  |    26.34867    12.77367      2.06    0.039      1.293414     51.40392
             |
       _cons |    67.66861    1.387883     48.76    0.000      64.94631     70.39091
-------------+----------------------------------------------------------------

. table vag_disch rp, statistic(mean wpc) nototal


---------------------------+--------------------------------
                           |   Retained placenta at calving
                           |              no             yes
---------------------------+--------------------------------
Vaginal discharge observed |
  no                       |        67.66861         74.0084
  yes                      |        68.21154           100.9
---------------------------+--------------------------------

. * table vag_disch rp, c(mean wpc) // version 16 or before
. testparm i.rp#i.vag_disch // same as t-test in coefficient table


 ( 1)  1.rp#1.vag_disch = 0


       F(  1,  1570) =    4.25
            Prob > F =    0.0393

. regress wpc i.rp##i.vag_disch // same with less typing

      Source |       SS           df       MS      Number of obs   =     1,574
-------------+----------------------------------   F(3, 1570)      =      4.53
       Model |  35915.9774         3  11971.9925   Prob > F        =    0.0036
    Residual |  4152174.58     1,570  2644.69719   R-squared       =    0.0086
-------------+----------------------------------   Adj R-squared   =    0.0067
       Total |  4188090.56     1,573  2662.48605   Root MSE        =    51.427

-------------+----------------------------------------------------------------
          wpc | Coefficient   Std. err.       t    P>|t|      [95% con f. interval]
-------------+----------------------------------------------------------------
          rp |
        yes  |    6.339794    4.914322      1.29    0.197     -3.299531     15.97912
             |
   vag_disch |
        yes  |    .5429296    7.265382      0.07    0.940     -13.70794      14.7938
             |
rp#vag_disch |
    yes#yes  |    26.34867    12.77367      2.06    0.039      1.293414     51.40392
             |
       _cons |    67.66861    1.387883     48.76    0.000      64.94631     70.39091
-------------+----------------------------------------------------------------

. regress wpc rp##vag_disch // i. factor notation may be omitted, but maybe keep it

      Source |       SS           df       MS      Number of obs   =     1,574
-------------+----------------------------------   F(3, 1570)      =      4.53
       Model |  35915.9774         3  11971.9925   Prob > F        =    0.0036
    Residual |  4152174.58     1,570  2644.69719   R-squared       =    0.0086
-------------+----------------------------------   Adj R-squared   =    0.0067
       Total |  4188090.56     1,573  2662.48605   Root MSE        =    51.427

-------------+----------------------------------------------------------------
          wpc | Coefficient   Std. err.       t    P>|t|      [95% con f. interval]
-------------+----------------------------------------------------------------
          rp |
        yes  |    6.339794    4.914322      1.29    0.197     -3.299531     15.97912
             |
   vag_disch |
        yes  |    .5429296    7.265382      0.07    0.940     -13.70794      14.7938
             |
rp#vag_disch |
    yes#yes  |    26.34867    12.77367      2.06    0.039      1.293414     51.40392
             |
       _cons |    67.66861    1.387883     48.76    0.000      64.94631     70.39091
-------------+----------------------------------------------------------------

. margins i.rp#i.vag_disch

 Adjusted predictions                                    Number of obs    = 1,574
 Model VCE: OLS


Expression:  Linear prediction, predict()


-------------+----------------------------------------------------------------
             |            Delta-method
             |     Margin    std. err.       t    P>|t|      [95% con f. interval]
-------------+----------------------------------------------------------------
rp#vag_disch |
      no#no  |    67.66861    1.387883     48.76    0.000      64.94631     70.39091
     no#yes  |    68.21154    7.131589      9.56    0.000       54.2231     82.19998
     yes#no  |     74.0084     4.71427     15.70    0.000      64.76147     83.25533
    yes#yes  |       100.9    9.389173     10.75    0.000      82.48336     119.3166
-------------+----------------------------------------------------------------

. marginsplot, noci /* interaction plot in Stata, more details in Lecture 3a */


Variables that uniquely identify margins: rp vag_disch


.
. * one continuous, one categorical predictor (VER 14.10)
. regress wpc i.dyst milk120

      Source |       SS           df       MS      Number of obs   =     1,536
-------------+----------------------------------   F(2, 1533)      =      1.01
       Model |   5380.7883         2  2690.39415   Prob > F        =    0.3658
    Residual |  4098983.87     1,533  2673.83162   R-squared       =    0.0013
-------------+----------------------------------   Adj R-squared   =    0.0000
       Total |  4104364.66     1,535   2673.8532   Root MSE        =    51.709

-------------+----------------------------------------------------------------
          wpc | Coefficient   Std. err.       t    P>|t|      [95% con f. interval]
-------------+----------------------------------------------------------------
        dyst |
        yes  |    4.004176    5.568429      0.72    0.472     -6.918368     14.92672
     milk120 |   -.0022375    .0018933     -1.18    0.237     -.0059512     .0014763
       _cons |    75.93882    6.255198     12.14    0.000      63.66917     88.20847
-------------+----------------------------------------------------------------

. generate milk120k= milk120/1000 // rescale predictor
(38 missing values generated)

. regress wpc i.dyst milk120k // additive model

      Source |       SS           df       MS      Number of obs   =     1,536
-------------+----------------------------------   F(2, 1533)      =      1.01
       Model |  5380.78796         2  2690.39398   Prob > F        =    0.3658
    Residual |  4098983.87     1,533  2673.83162   R-squared       =    0.0013
-------------+----------------------------------   Adj R-squared   =    0.0000
       Total |  4104364.66     1,535   2673.8532   Root MSE        =    51.709

-------------+----------------------------------------------------------------
          wpc | Coefficient   Std. err.       t    P>|t|      [95% con f. interval]
-------------+----------------------------------------------------------------
        dyst |
        yes  |    4.004176    5.568429      0.72    0.472     -6.918368     14.92672
    milk120k |   -2.237454    1.893316     -1.18    0.237     -5.951217     1.476309
       _cons |    75.93882    6.255198     12.14    0.000      63.66917     88.20847
-------------+----------------------------------------------------------------

. regress wpc i.dyst milk120k i.dyst#c.milk120k // interaction added, see 2aL for factor notation

      Source |       SS           df       MS      Number of obs   =     1,536
-------------+----------------------------------   F(3, 1532)      =      3.83
       Model |   30572.875         3  10190.9583   Prob > F        =    0.0095
    Residual |  4073791.78     1,532  2659.13302   R-squared       =    0.0074
-------------+----------------------------------   Adj R-squared   =    0.0055
       Total |  4104364.66     1,535   2673.8532   Root MSE        =    51.567

----------------+----------------------------------------------------------------
             wpc | Coefficient   Std. err.       t    P>|t|      [95% con f. interval]
----------------+----------------------------------------------------------------
           dyst |
           yes  |   -85.48838    29.60089     -2.89    0.004     -143.5509    -27.42583
       milk120k |   -3.446531    1.928535     -1.79    0.074     -7.229379     .3363161
                |
dyst#c.milk120k |
           yes  |    29.14238    9.468101      3.08    0.002      10.57056     47.71419
                |
          _cons |    79.83774    6.365297     12.54    0.000      67.35213     92.32336
----------------+----------------------------------------------------------------

. regress wpc i.dyst##c.milk120k

      Source |       SS           df       MS      Number of obs   =     1,536
-------------+----------------------------------   F(3, 1532)      =      3.83
       Model |   30572.875         3  10190.9583   Prob > F        =    0.0095
    Residual |  4073791.78     1,532  2659.13302   R-squared       =    0.0074
-------------+----------------------------------   Adj R-squared   =    0.0055
       Total |  4104364.66     1,535   2673.8532   Root MSE        =    51.567

----------------+----------------------------------------------------------------
             wpc | Coefficient   Std. err.       t    P>|t|      [95% con f. interval]
----------------+----------------------------------------------------------------
           dyst |
           yes  |   -85.48838    29.60089     -2.89    0.004     -143.5509    -27.42583
       milk120k |   -3.446531    1.928535     -1.79    0.074     -7.229379     .3363161
                |
dyst#c.milk120k |
           yes  |    29.14238    9.468101      3.08    0.002      10.57056     47.71419
                |
          _cons |    79.83774    6.365297     12.54    0.000      67.35213     92.32336
----------------+----------------------------------------------------------------

. margins i.dyst, at(milk120k=(1(1)6))

 Adjusted predictions                                    Number of obs    = 1,536
 Model VCE: OLS


Expression:  Linear prediction, predict()


1._at:   milk120k   = 1
2._at:   milk120k   = 2
3._at:   milk120k   = 3
4._at:   milk120k   = 4
5._at:   milk120k   = 5
6._at:   milk120k   = 6

-------------+----------------------------------------------------------------
             |            Delta-method
             |     Margin    std. err.       t    P>|t|      [95% con f. interval]
-------------+----------------------------------------------------------------
    _at#dyst |
       1#no  |    76.39121    4.499919     16.98    0.000      67.56456     85.21787
      1#yes  |    20.04521     19.8754      1.01    0.313     -18.94067     59.03108
       2#no  |    72.94468    2.723976     26.78    0.000      67.60156      78.2878
      2#yes  |    45.74105    11.23473      4.07    0.000      23.70398     67.77812
       3#no  |    69.49815    1.424536     48.79    0.000       66.7039      72.2924
      3#yes  |    71.43689    5.409061     13.21    0.000      60.82695     82.04684
       4#no  |    66.05162    2.019173     32.71    0.000      62.09098     70.01225
      4#yes  |    97.13274     10.2053      9.52    0.000       77.1149     117.1506
       5#no  |    62.60509    3.682839     17.00    0.000      55.38115     69.82903
      5#yes  |    122.8286    18.73205      6.56    0.000       86.0854     159.5718
       6#no  |    59.15855    5.521598     10.71    0.000      48.32786     69.98925
      6#yes  |    148.5244    27.73955      5.35    0.000      94.11291     202.9359
-------------+----------------------------------------------------------------

. marginsplot, noci /* interaction plot in Stata, more details in Lecture 3a */


Variables that uniquely identify margins: milk120k dyst


.
. * two continuous predictors (VER 14.11)
. regress wpc c.parity##c.milk120k

      Source |       SS           df       MS      Number of obs   =     1,536
-------------+----------------------------------   F(3, 1532)      =      2.26
       Model |  18084.1905         3  6028.06351   Prob > F        =    0.0797
    Residual |  4086280.47     1,532   2667.2849   R-squared       =    0.0044
-------------+----------------------------------   Adj R-squared   =    0.0025
       Total |  4104364.66     1,535   2673.8532   Root MSE        =    51.646

--------------------+----------------------------------------------------------------
                 wpc | Coefficient   Std. err.       t    P>|t|      [95% con f. interval]
--------------------+----------------------------------------------------------------
             parity |    4.889989    4.437817      1.10    0.271     -3.814849     13.59483
           milk120k |   -1.666398    4.212531     -0.40    0.692     -9.929335     6.596539
                    |
c.parity#c.milk120k |   -.8764359    1.363504     -0.64    0.520     -3.550968     1.798096
                    |
              _cons |    69.02201    12.77659      5.40    0.000      43.96056     94.08347
--------------------+----------------------------------------------------------------

. * interaction plot in Lecture 3a
.
. * illustration of confounding (not a VER example)
. regress wpc i.vag_disch // crude model: vag_disch adds ~12 days, P=0.04

      Source |       SS           df       MS      Number of obs   =     1,574
-------------+----------------------------------   F(1, 1572)      =      4.21
       Model |  11186.2576         1  11186.2576   Prob > F        =    0.0404
    Residual |   4176904.3     1,572   2657.0638   R-squared       =    0.0027
-------------+----------------------------------   Adj R-squared   =    0.0020
       Total |  4188090.56     1,573  2662.48605   Root MSE        =    51.547

-------------+----------------------------------------------------------------
          wpc | Coefficient   Std. err.       t    P>|t|      [95% con f. interval]
-------------+----------------------------------------------------------------
   vag_disch |
        yes  |    11.99647    5.846716      2.05    0.040      .5282858     23.46465
       _cons |    68.17426    1.334494     51.09    0.000      65.55669     70.79184
-------------+----------------------------------------------------------------

. regress wpc i.herd if vag_disch==0 // confounder and outcome in unexposed individuals

      Source |       SS           df       MS      Number of obs   =     1,492
-------------+----------------------------------   F(6, 1485)      =     14.92
       Model |  223706.812         6  37284.4686   Prob > F        =    0.0000
    Residual |  3711781.88     1,485  2499.51642   R-squared       =    0.0568
-------------+----------------------------------   Adj R-squared   =    0.0530
       Total |  3935488.69     1,491   2639.4961   Root MSE        =    49.995

-------------+----------------------------------------------------------------
          wpc | Coefficient   Std. err.       t    P>|t|      [95% con f. interval]
-------------+----------------------------------------------------------------
        herd |
          2  |   -7.455598    4.572715     -1.63    0.103     -16.42527      1.51407
          3  |    12.30396    4.184586      2.94    0.003      4.095636     20.51229
          4  |    -20.0733     4.70619     -4.27    0.000     -29.30479    -10.84181
          5  |   -21.78125    5.489853     -3.97    0.000     -32.54994    -11.01256
        106  |   -15.40129    4.618509     -3.33    0.001     -24.46079    -6.341796
        119  |   -17.26021     5.05679     -3.41    0.001     -27.17942       -7.341
             |
       _cons |     75.1583    3.106548     24.19    0.000      69.06461     81.25199
-------------+----------------------------------------------------------------

. testparm i.herd // significant herd effects


 ( 1)  2.herd = 0


 ( 2)  3.herd = 0


 ( 3)  4.herd = 0


 ( 4)  5.herd = 0


 ( 5)  106.herd = 0


 ( 6)  119.herd = 0


       F(  6,  1485) =   14.92
            Prob > F =    0.0000

. tabulate herd vag_disch, chi2 row // vag_disch varies across significantly herds

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

           |   Vaginal discharge
      Herd |       observed
    Number |        no        yes |     Total
-----------+----------------------+----------
         1 |       259         13 |       272
           |     95.22       4.78 |    100.00
-----------+----------------------+----------
         2 |       222          5 |       227
           |     97.80       2.20 |    100.00
-----------+----------------------+----------
         3 |       318          4 |       322
           |     98.76       1.24 |    100.00
-----------+----------------------+----------
         4 |       200          3 |       203
           |     98.52       1.48 |    100.00
-----------+----------------------+----------
         5 |       122          7 |       129
           |     94.57       5.43 |    100.00
-----------+----------------------+----------
       106 |       214         39 |       253
           |     84.58      15.42 |    100.00
-----------+----------------------+----------
       119 |       157         11 |       168
           |     93.45       6.55 |    100.00
-----------+----------------------+----------
     Total |     1,492         82 |     1,574
           |     94.79       5.21 |    100.00

          Pearson chi2(6) =  74.2267   Pr = 0.000

. regress wpc i.vag_disch i.herd // adjusted model

      Source |       SS           df       MS      Number of obs   =     1,574
-------------+----------------------------------   F(7, 1566)      =     14.35
       Model |   252509.59         7  36072.7985   Prob > F        =    0.0000
    Residual |  3935580.97     1,566  2513.14238   R-squared       =    0.0603
-------------+----------------------------------   Adj R-squared   =    0.0561
       Total |  4188090.56     1,573  2662.48605   Root MSE        =    50.131

-------------+----------------------------------------------------------------
          wpc | Coefficient   Std. err.       t    P>|t|      [95% con f. interval]
-------------+----------------------------------------------------------------
   vag_disch |
        yes  |    17.81936    5.825177      3.06    0.002      6.393393     29.24533
             |
        herd |
          2  |   -8.178615    4.509228     -1.81    0.070     -17.02338     .6661455
          3  |    11.71303    4.133611      2.83    0.005      3.605039     19.82103
          4  |   -20.74627    4.653654     -4.46    0.000     -29.87432    -11.61823
          5  |   -22.15954    5.359351     -4.13    0.000      -32.6718    -11.64728
        106  |   -18.18881    4.422295     -4.11    0.000     -26.86305    -9.514562
        119  |   -17.85657    4.920293     -3.63    0.000     -27.50763    -8.205518
             |
       _cons |    75.97555    3.052377     24.89    0.000      69.98837     81.96272
-------------+----------------------------------------------------------------

. display "relative change =" abs(12.00-17.82)/12.00
relative change =.485

.
. * VER 14.12 multivariable causal model
. use daisy2red, clear

. generate calv_mth=month(calv_dt)

. tabulate calv_mth, summ(wpc)

            |    Summary of Interval from wait
            |        period to conception
   calv_mth |        Mean   Std. dev.       Freq.
------------+------------------------------------
          1 |    77.00885   50.850533         113
          2 |   60.009615   45.765017         104
          3 |   69.402516   54.240687         159
          4 |   60.385827   53.046258         127
          5 |          63   49.882911         106
          6 |   64.350365   51.665154         137
          7 |   63.581818   45.047453         110
          8 |   69.514019   54.277478         107
          9 |   74.664773   52.827873         176
         10 |   72.477011   49.332856         174
         11 |   72.026846   54.231479         149
         12 |   73.508929   53.222801         112
------------+------------------------------------
      Total |   68.799238   51.599283       1,574

. generate aut_calv=(calv_mth>=2 & calv_mth<=7) if !missing(calv_mth)

. tabulate aut_calv, summ(wpc)

            |    Summary of Interval from wait
            |        period to conception
   aut_calv |        Mean   Std. dev.       Freq.
------------+------------------------------------
          0 |   73.233454   52.236588         831
          1 |   63.839838   50.451973         743
------------+------------------------------------
      Total |   68.799238   51.599283       1,574

. generate hs100=herd_size/100

. generate parity1=parity-1

. regress wpc c.hs100##c.hs100 parity1 i.aut_calv i.twin i.dyst i.rp##i.vag_disch

      Source |       SS           df       MS      Number of obs   =     1,574
-------------+----------------------------------   F(9, 1564)      =     13.22
       Model |  296062.681         9  32895.8535   Prob > F        =    0.0000
    Residual |  3892027.88     1,564  2488.50887   R-squared       =    0.0707
-------------+----------------------------------   Adj R-squared   =    0.0653
       Total |  4188090.56     1,573  2662.48605   Root MSE        =    49.885

----------------+----------------------------------------------------------------
             wpc | Coefficient   Std. err.       t    P>|t|      [95% con f. interval]
----------------+----------------------------------------------------------------
          hs100 |   -36.05705    15.05032     -2.40    0.017     -65.57798     -6.53612
                |
c.hs100#c.hs100 |    11.13827    3.111145      3.58    0.000      5.035818     17.24073
                |
        parity1 |     1.13721    .8583103      1.32    0.185       -.54635      2.82077
     1.aut_calv |   -8.263839    2.537751     -3.26    0.001     -13.24159    -3.286086
                |
           twin |
           yes  |    20.68314    9.845165      2.10    0.036       1.37203     39.99425
                |
           dyst |
           yes  |    11.70041    5.462576      2.14    0.032      .9856659     22.41516
                |
             rp |
           yes  |     5.98687    4.811976      1.24    0.214     -3.451734     15.42547
                |
      vag_disch |
           yes  |    1.228195    7.161395      0.17    0.864     -12.81875     15.27514
                |
   rp#vag_disch |
       yes#yes  |    22.85194    12.51605      1.83    0.068     -1.698056     47.40194
                |
          _cons |    84.66125    17.61671      4.81    0.000      50.10639     119.2161
----------------+----------------------------------------------------------------

. * same results as in VER results, except for intercept
. * to center herdsize at 251 is pretty meaningless...
. summarize herd_size, d // this is not the actual distribution of herd_size!

                          Herd size
-------------------------------------------------------------
      Percentiles      Smallest
 1%          125            125
 5%          125            125
10%          185            125       Obs               1,574
25%          201            125       Sum of wgt.       1,574

50%          263                      Mean           251.0076
                        Largest       Std. dev.      62.01692
75%          294            333
90%          333            333       Variance       3846.098
95%          333            333       Skewness      -.3550929
99%          333            333       Kurtosis       2.256969

. tabulate herd_size

  Herd size |      Freq.     Percent        Cum.
------------+-----------------------------------
        125 |        129        8.20        8.20
        185 |        168       10.67       18.87
        201 |        203       12.90       31.77
        235 |        227       14.42       46.19
        263 |        253       16.07       62.26
        294 |        272       17.28       79.54
        333 |        322       20.46      100.00
------------+-----------------------------------
      Total |      1,574      100.00

. bysort herd: gen within=_n // generates counter 1,2,... of cows in each herd

. summarize herd_size if within==1,d

                          Herd size
-------------------------------------------------------------
      Percentiles      Smallest
 1%          125            125
 5%          125            185
10%          125            201       Obs                   7
25%          185            235       Sum of wgt.           7

50%          235                      Mean           233.7143
                        Largest       Std. dev.      70.25362
75%          294            235
90%          333            263       Variance       4935.571
95%          333            294       Skewness      -.1103475
99%          333            333       Kurtosis       2.056033

.
end of do-file

. exit, clear