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

. * do-file for lecture 3a of VHM 802/812, Winter 2025
. version 18 /* works also with versions 14-17, except for table command */

. set more off

. set scheme stcolor_alt

. cd "r:\"
r:\

.
. use micegrp.dta, clear

. regress dead dose /* to illustrate extrapolation below zero in regression */

      Source |       SS           df       MS      Number of obs   =        12
-------------+----------------------------------   F(1, 10)        =     29.83
       Model |   49.616367         1   49.616367   Prob > F        =    0.0003
    Residual |   16.633633        10   1.6633633   R-squared       =    0.7489
-------------+----------------------------------   Adj R-squared   =    0.7238
       Total |       66.25        11  6.02272727   Root MSE        =    1.2897

-------------+----------------------------------------------------------------
         dead | Coefficient   Std. err.       t    P>|t|      [95% con f. interval]
-------------+----------------------------------------------------------------
        dose |     29.1441    5.336195      5.46    0.000      17.25431     41.03388
       _cons |   -2.107171    1.049163     -2.01    0.072     -4.444851     .2305092
-------------+----------------------------------------------------------------

. use mice.dta, clear

. regress dead dose /* same estimates of regression coefficients, up to scaling */

      Source |       SS           df       MS      Number of obs   =       120
-------------+----------------------------------   F(1, 118)       =     27.41
       Model |   4.9616367         1   4.9616367   Prob > F        =    0.0000
    Residual |  21.3633633       118  .181045452   R-squared       =    0.1885
-------------+----------------------------------   Adj R-squared   =    0.1816
       Total |      26.325       119  .221218487   Root MSE        =    .42549

-------------+----------------------------------------------------------------
         dead | Coefficient   Std. err.       t    P>|t|      [95% con f. interval]
-------------+----------------------------------------------------------------
        dose |     2.91441    .5567137      5.24    0.000      1.811965     4.016855
       _cons |   -.2107171    .1094569     -1.93    0.057     -.4274715     .0060373
-------------+----------------------------------------------------------------

.
. * logistic regression using logit command
. logit dead dose

Iteration 0:  Log likelihood = -75.669723
Iteration 1:  Log likelihood = -64.201632
Iteration 2:  Log likelihood = -63.944944
Iteration 3:  Log likelihood = -63.944713
Iteration 4:  Log likelihood = -63.944713

 Logistic regression                                    Number of obs    = 120
                                                        LR chi2(1)    = 23.45
                                                        Prob > chi2    = 0.0000
 Log likelihood    = -63.944713                         Pseudo R2     = 0.1549

-------------+----------------------------------------------------------------
         dead | Coefficient   Std. err.       z    P>|z|      [95% con f. interval]
-------------+----------------------------------------------------------------
        dose |     14.6369    3.332533      4.39    0.000       8.10526     21.16855
       _cons |   -3.569738    .7053432     -5.06    0.000     -4.952185    -2.187291
-------------+----------------------------------------------------------------

. predict linpred, xb /* predictions on logit scale */

. twoway (scatter linpred dose) (line linpred dose)

. lincom _cons+dose*0.1


 ( 1)  .1*[dead]dose + [dead]_cons = 0


-------------+----------------------------------------------------------------
         dead | Coefficient   Std. err.       z    P>|z|      [95% con f. interval]
-------------+----------------------------------------------------------------
         (1) |   -2.106048    .4015256     -5.25    0.000     -2.893023    -1.319072
-------------+----------------------------------------------------------------

. predict phat, p     /* predicted probabilities */

. twoway (scatter phat dose) (line phat dose)

. estat gof           /* Pearson goodness-of-fit test */

Goodness-of-fit test after logistic model
Variable: dead

Number of observations   = 120
Number of covariate patterns   = 12
Pearson chi2(10)   = 8.74
Prob > chi2   = 0.5567

.
. * analysis with dichotomous dose
. prtesti 60 32 60 7, count /* 2 independent proportions */

Two-sample test of proportions                     x: Number of obs =       60
                                                   y: Number of obs =       60
-------------+----------------------------------------------------------------
             |       Mean   Std. err.      z    P>|z|     [95% conf. interval]
-------------+----------------------------------------------------------------
           x |   .5333333   .0644061                      .4070997     .659567
           y |   .1166667   .0414438                      .0354382    .1978951
-------------+----------------------------------------------------------------
        diff |   .4166667   .0765881                      .2665567    .5667766
             |  under H0:   .0855132     4.87   0.000
-------------+----------------------------------------------------------------
        diff = prop(x) - prop(y)                                  z =   4.8725
    H0: diff = 0

    Ha: diff < 0                 Ha: diff != 0                 Ha: diff > 0
 Pr(Z < z) = 1.0000         Pr(|Z| > |z|) = 0.0000          Pr(Z > z) = 0.0000

. csi 32 7 28 53, or        /* relative risk and odds-ratio */

                 |   Exposed   Unexposed  |      Total
-----------------+------------------------+-----------
           Cases |        32           7  |         39
        Noncases |        28          53  |         81
-----------------+------------------------+-----------
           Total |        60          60  |        120
                 |                        |
            Risk |  .5333333    .1166667  |       .325
                 |                        |
                 |      Point estimate    |    [95% conf. interval]
                 +------------------------+------------------------
 Risk difference |         .4166667       |    .2665567    .5667766
      Risk ratio |         4.571429       |    2.191203     9.53721
 Attr. frac. ex. |           .78125       |    .5436296    .8951475
 Attr. frac. pop |         .6410256       |
      Odds ratio |         8.653061       |    3.441938    21.63764 (Cornfield)
                 +------------------------+------------------------
                               chi2(1) =    23.74  Pr>chi2 = 0.0000

. generate dose2=(dose>0.16) if dose~=.

. tabulate dead dose2, chi2 column exact lrchi2 /* 2x2 table analysis */

+-------------------+
| Key               |
+-------------------+
|     frequency     |
| column percentage |
+-------------------+

           |         dose2
      dead |         0          1 |     Total
-----------+----------------------+----------
         0 |        53         28 |        81
           |     88.33      46.67 |     67.50
-----------+----------------------+----------
         1 |         7         32 |        39
           |     11.67      53.33 |     32.50
-----------+----------------------+----------
     Total |        60         60 |       120
           |    100.00     100.00 |    100.00

          Pearson chi2(1) =  23.7417   Pr = 0.000
 Likelihood-ratio chi2(1) =  25.2010   Pr = 0.000
           Fisher's exact =                 0.000
   1-sided Fisher's exact =                 0.000

. logit dead i.dose2    /* logistic regression */

Iteration 0:  Log likelihood = -75.669723
Iteration 1:  Log likelihood = -63.522292
Iteration 2:  Log likelihood = -63.070143
Iteration 3:  Log likelihood =  -63.06923
Iteration 4:  Log likelihood =  -63.06923

 Logistic regression                                    Number of obs    = 120
                                                        LR chi2(1)    = 25.20
                                                        Prob > chi2    = 0.0000
 Log likelihood    = -63.06923                          Pseudo R2     = 0.1665

-------------+----------------------------------------------------------------
         dead | Coefficient   Std. err.       z    P>|z|      [95% con f. interval]
-------------+----------------------------------------------------------------
     1.dose2 |    2.157913    .4782148      4.51    0.000      1.220629     3.095197
       _cons |   -2.024382    .4021506     -5.03    0.000     -2.812582    -1.236181
-------------+----------------------------------------------------------------

. logistic dead i.dose2 /* same as logit with option or */

 Logistic regression                                    Number of obs    = 120
                                                        LR chi2(1)    = 25.20
                                                        Prob > chi2    = 0.0000
 Log likelihood    = -63.06923                          Pseudo R2     = 0.1665

-------------+----------------------------------------------------------------
         dead | Odds ratio    Std. err.       z    P>|z|      [95% con f. interval]
-------------+----------------------------------------------------------------
     1.dose2 |    8.653061    4.138022      4.51    0.000       3.38932     22.09159
       _cons |    .1320755    .0531142     -5.03    0.000      .0600497     .2904914
-------------+----------------------------------------------------------------


Note: _cons estimates baseline odds.


. display 1/(1+exp(2.024382)) /* baseline prob */
.11666664

. display invlogit(-2.024382) /* same with invlogit function */
.11666664

. logit dead dose

Iteration 0:  Log likelihood = -75.669723
Iteration 1:  Log likelihood = -64.201632
Iteration 2:  Log likelihood = -63.944944
Iteration 3:  Log likelihood = -63.944713
Iteration 4:  Log likelihood = -63.944713

 Logistic regression                                    Number of obs    = 120
                                                        LR chi2(1)    = 23.45
                                                        Prob > chi2    = 0.0000
 Log likelihood    = -63.944713                         Pseudo R2     = 0.1549

-------------+----------------------------------------------------------------
         dead | Coefficient   Std. err.       z    P>|z|      [95% con f. interval]
-------------+----------------------------------------------------------------
        dose |     14.6369    3.332533      4.39    0.000       8.10526     21.16855
       _cons |   -3.569738    .7053432     -5.06    0.000     -4.952185    -2.187291
-------------+----------------------------------------------------------------

. logistic dead dose  /* same as logit with option or */

 Logistic regression                                    Number of obs    = 120
                                                        LR chi2(1)    = 23.45
                                                        Prob > chi2    = 0.0000
 Log likelihood    = -63.944713                         Pseudo R2     = 0.1549

-------------+----------------------------------------------------------------
         dead | Odds ratio    Std. err.       z    P>|z|      [95% con f. interval]
-------------+----------------------------------------------------------------
        dose |     2273668     7577074      4.39    0.000      3311.843     1.56e+09
       _cons |    .0281632    .0198647     -5.06    0.000      .0070679     .1122204
-------------+----------------------------------------------------------------


Note: _cons estimates baseline odds.


. display 1/(1+exp(3.569738)) /* intercept prob */
.02739179

. display invlogit(-3.569738) /* same with invlogit function */
.02739179

.
. * demonstration of analysis using grouped data
. use micegrp.dta, clear

. glm dead dose, fam(binomial total)

Iteration 0:  Log likelihood =  -18.60686
Iteration 1:  Log likelihood = -18.516941
Iteration 2:  Log likelihood = -18.516789
Iteration 3:  Log likelihood = -18.516789

Generalized linear models                         Number of obs   =         12
Optimization     : ML                             Residual df     =         10
                                                  Scale parameter =          1
Deviance         =  10.25371404                   (1/df) Deviance =   1.025371
Pearson          =    8.7421447                   (1/df) Pearson  =   .8742145

Variance function: V(u) =  u*(1-u/total)          [Binomial]
Link function    : g(u) =  ln(u/(total-u))        [Logit]

                                                  AIC             =   3.419465
Log likelihood   = -18.51678904                   BIC             =  -14.59535

-------------+----------------------------------------------------------------
             |                 OIM
         dead | Coefficient   std. err.       z    P>|z|      [95% con f. interval]
-------------+----------------------------------------------------------------
        dose |     14.6369    3.332533      4.39    0.000       8.10526     21.16855
       _cons |   -3.569738    .7053432     -5.06    0.000     -4.952185    -2.187291
-------------+----------------------------------------------------------------

. predict linpred, xb

. predict deadhat, mu

. generate phat=deadhat/total /* or as: invlogit(linpred) */

. * Pearson goodness-of-fit test
. display chi2tail(12-1-1,8.7421447)
.55673742

. display chi2tail(12-1-1,e(deviance_p))
.55673742

.
. * prediction part of lecture
. use daisy2red, clear

.
. * 1a)
. regress wpc i.herd

      Source |       SS           df       MS      Number of obs   =     1,574
-------------+----------------------------------   F(6, 1567)      =     15.11
       Model |  228992.526         6   38165.421   Prob > F        =    0.0000
    Residual |  3959098.03     1,567  2526.54629   R-squared       =    0.0547
-------------+----------------------------------   Adj R-squared   =    0.0511
       Total |  4188090.56     1,573  2662.48605   Root MSE        =    50.265

-------------+----------------------------------------------------------------
          wpc | Coefficient   Std. err.       t    P>|t|      [95% con f. interval]
-------------+----------------------------------------------------------------
        herd |
          2  |   -8.637779    4.518731     -1.91    0.056     -17.50118     .2256182
          3  |    11.08273    4.139467      2.68    0.007      2.963254     19.20221
          4  |    -21.3346    4.662061     -4.58    0.000     -30.47913    -12.19006
          5  |   -22.04426    5.373491     -4.10    0.000     -32.58425    -11.50427
        106  |   -16.29361    4.390344     -3.71    0.000     -24.90518    -7.682041
        119  |   -17.54149    4.932316     -3.56    0.000     -27.21613    -7.866858
             |
       _cons |    76.82721    3.047749     25.21    0.000      70.84911      82.8053
-------------+----------------------------------------------------------------

. margins herd

 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]
-------------+----------------------------------------------------------------
        herd |
          1  |    76.82721    3.047749     25.21    0.000      70.84911      82.8053
          2  |    68.18943    3.336189     20.44    0.000      61.64556     74.73329
          3  |    87.90994    2.801146     31.38    0.000      82.41555     93.40433
          4  |    55.49261    3.527895     15.73    0.000      48.57272      62.4125
          5  |    54.78295    4.425566     12.38    0.000      46.10229      63.4636
        106  |     60.5336    3.160119     19.16    0.000      54.33509      66.7321
        119  |    59.28571    3.878011     15.29    0.000      51.67908     66.89235
-------------+----------------------------------------------------------------

. lincom _cons+2.herd


 ( 1)  2.herd + _cons = 0


-------------+----------------------------------------------------------------
          wpc | Coefficient   Std. err.       t    P>|t|      [95% con f. interval]
-------------+----------------------------------------------------------------
         (1) |    68.18943    3.336189     20.44    0.000      61.64556     74.73329
-------------+----------------------------------------------------------------

. tabstat wpc, statistics(mean sd semean count) by(herd)

Summary for variables: wpc
Group variable: herd (Herd Number)

herd   |      Mean        SD  se(mean)         N
---------+----------------------------------------
1   |  76.82721  54.25425  3.289647       272
2   |  68.18943  51.42065  3.412909       227
3   |  87.90994  56.03043  3.122454       322
4   |  55.49261  46.33257  3.251909       203
5   |  54.78295  43.38904  3.820192       129
106   |   60.5336  43.14395  2.712437       253
119   |  59.28571  49.69417  3.833989       168
---------+----------------------------------------
Total   |  68.79924  51.59928  1.300593      1574
---------+----------------------------------------

.
. * 1b)
. regress wpc milk120

      Source |       SS           df       MS      Number of obs   =     1,536
-------------+----------------------------------   F(1, 1534)      =      1.50
       Model |  3998.19335         1  3998.19335   Prob > F        =    0.2215
    Residual |  4100366.46     1,534  2672.98987   R-squared       =    0.0010
-------------+----------------------------------   Adj R-squared   =    0.0003
       Total |  4104364.66     1,535   2673.8532   Root MSE        =    51.701

-------------+----------------------------------------------------------------
          wpc | Coefficient   Std. err.       t    P>|t|      [95% con f. interval]
-------------+----------------------------------------------------------------
     milk120 |   -.0023117    .0018902     -1.22    0.222     -.0060194     .0013959
       _cons |    76.41752    6.218695     12.29    0.000      64.21947     88.61556
-------------+----------------------------------------------------------------

. margins , at( milk120=(1400(1000)5400) )

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


Expression:  Linear prediction, predict()


1._at:   milk120   = 1400
2._at:   milk120   = 2400
3._at:   milk120   = 3400
4._at:   milk120   = 4400
5._at:   milk120   = 5400

-------------+----------------------------------------------------------------
             |            Delta-method
             |     Margin    std. err.       t    P>|t|      [95% con f. interval]
-------------+----------------------------------------------------------------
         _at |
          1  |    73.18107    3.675762     19.91    0.000      65.97102     80.39112
          2  |    70.86932    2.028289     34.94    0.000      66.89081     74.84783
          3  |    68.55757    1.364691     50.24    0.000      65.88072     71.23443
          4  |    66.24583    2.599327     25.49    0.000      61.14722     71.34444
          5  |    63.93408    4.335471     14.75    0.000         55.43     72.43816
-------------+----------------------------------------------------------------

. marginsplot


Variables that uniquely identify margins: milk120


. margins , at( milk120=(1400(1000)5400) ) plot /* same as with marginsplot, but less options */

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


Expression:  Linear prediction, predict()


1._at:   milk120   = 1400
2._at:   milk120   = 2400
3._at:   milk120   = 3400
4._at:   milk120   = 4400
5._at:   milk120   = 5400

-------------+----------------------------------------------------------------
             |            Delta-method
             |     Margin    std. err.       t    P>|t|      [95% con f. interval]
-------------+----------------------------------------------------------------
         _at |
          1  |    73.18107    3.675762     19.91    0.000      65.97102     80.39112
          2  |    70.86932    2.028289     34.94    0.000      66.89081     74.84783
          3  |    68.55757    1.364691     50.24    0.000      65.88072     71.23443
          4  |    66.24583    2.599327     25.49    0.000      61.14722     71.34444
          5  |    63.93408    4.335471     14.75    0.000         55.43     72.43816
-------------+----------------------------------------------------------------


Variables that uniquely identify margins: milk120


. lincom _cons+milk120*3400


 ( 1)  3400*milk120 + _cons = 0


-------------+----------------------------------------------------------------
          wpc | Coefficient   Std. err.       t    P>|t|      [95% con f. interval]
-------------+----------------------------------------------------------------
         (1) |    68.55757    1.364691     50.24    0.000      65.88072     71.23443
-------------+----------------------------------------------------------------

. sum milk120, d

       Milk volume (l) in first 120 days of lactation
-------------------------------------------------------------
      Percentiles      Smallest
 1%       1698.5         1110.2
 5%       2077.2         1397.7
10%       2298.2         1461.5       Obs               1,536
25%      2731.95           1467       Sum of wgt.       1,536

50%      3215.25                      Mean           3215.096
                        Largest       Std. dev.      698.1316
75%       3682.1         5181.8
90%       4080.7           5278       Variance       487387.7
95%       4403.3         5399.9       Skewness       .1101838
99%       4904.4         5630.3       Kurtosis       2.845637

.
. * 1c)
. regress wpc c.milk120##c.milk120

      Source |       SS           df       MS      Number of obs   =     1,536
-------------+----------------------------------   F(2, 1533)      =      2.84
       Model |   15168.437         2  7584.21849   Prob > F        =    0.0585
    Residual |  4089196.22     1,533  2667.44698   R-squared       =    0.0037
-------------+----------------------------------   Adj R-squared   =    0.0024
       Total |  4104364.66     1,535   2673.8532   Root MSE        =    51.647

--------------------+----------------------------------------------------------------
                 wpc | Coefficient   Std. err.       t    P>|t|      [95% con f. interval]
--------------------+----------------------------------------------------------------
            milk120 |   -.0289184    .0131383     -2.20    0.028     -.0546894    -.0031475
                    |
c.milk120#c.milk120 |    4.09e-06    2.00e-06      2.05    0.041      1.70e-07     8.01e-06
                    |
              _cons |    117.7029    21.10975      5.58    0.000      76.29586     159.1099
--------------------+----------------------------------------------------------------

. margins , at( milk120=(1400(1000)5400) ) plot

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


Expression:  Linear prediction, predict()


1._at:   milk120   = 1400
2._at:   milk120   = 2400
3._at:   milk120   = 3400
4._at:   milk120   = 4400
5._at:   milk120   = 5400

-------------+----------------------------------------------------------------
             |            Delta-method
             |     Margin    std. err.       t    P>|t|      [95% con f. interval]
-------------+----------------------------------------------------------------
         _at |
          1  |     85.2313    6.939656     12.28    0.000      71.61908     98.84353
          2  |    71.85061    2.082156     34.51    0.000      67.76643     75.93478
          3  |    66.64766    1.652151     40.34    0.000      63.40695     69.88838
          4  |    69.62247    3.076558     22.63    0.000      63.58776     75.65717
          5  |    80.77502     9.29973      8.69    0.000      62.53348     99.01656
-------------+----------------------------------------------------------------


Variables that uniquely identify margins: milk120


. lincom _cons+milk120*3400+c.milk120#c.milk120*11560000


 ( 1)  3400*milk120 + 1.16e+07*c.milk120#c.milk120 + _cons = 0


-------------+----------------------------------------------------------------
          wpc | Coefficient   Std. err.       t    P>|t|      [95% con f. interval]
-------------+----------------------------------------------------------------
         (1) |    66.64766    1.652151     40.34    0.000      63.40695     69.88838
-------------+----------------------------------------------------------------

.
. * 1d)
. generate lnwpc=ln(wpc)

. regress lnwpc milk120

      Source |       SS           df       MS      Number of obs   =     1,536
-------------+----------------------------------   F(1, 1534)      =      2.51
       Model |  1.45963987         1  1.45963987   Prob > F        =    0.1130
    Residual |  890.377121     1,534  .580428371   R-squared       =    0.0016
-------------+----------------------------------   Adj R-squared   =    0.0010
       Total |  891.836761     1,535  .581001147   Root MSE        =    .76186

-------------+----------------------------------------------------------------
        lnwpc | Coefficient   Std. err.       t    P>|t|      [95% con f. interval]
-------------+----------------------------------------------------------------
     milk120 |   -.0000442    .0000279     -1.59    0.113     -.0000988     .0000105
       _cons |    4.103887    .0916379     44.78    0.000      3.924139     4.283636
-------------+----------------------------------------------------------------

. margins , at( milk120=(1400(1000)5400))

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


Expression:  Linear prediction, predict()


1._at:   milk120   = 1400
2._at:   milk120   = 2400
3._at:   milk120   = 3400
4._at:   milk120   = 4400
5._at:   milk120   = 5400

-------------+----------------------------------------------------------------
             |            Delta-method
             |     Margin    std. err.       t    P>|t|      [95% con f. interval]
-------------+----------------------------------------------------------------
         _at |
          1  |    4.042049    .0541655     74.62    0.000      3.935802     4.148295
          2  |    3.997878    .0298886    133.76    0.000      3.939252     4.056505
          3  |    3.953708    .0201099    196.60    0.000      3.914262     3.993154
          4  |    3.909538    .0383033    102.07    0.000      3.834405      3.98467
          5  |    3.865367    .0638869     60.50    0.000      3.740052     3.990682
-------------+----------------------------------------------------------------

. margins , at( milk120=(1400(1000)5400)) expression(exp(predict(xb))) plot

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


Expression:  exp(predict(xb))


1._at:   milk120   = 1400
2._at:   milk120   = 2400
3._at:   milk120   = 3400
4._at:   milk120   = 4400
5._at:   milk120   = 5400

-------------+----------------------------------------------------------------
             |            Delta-method
             |     Margin    std. err.       z    P>|z|      [95% con f. interval]
-------------+----------------------------------------------------------------
         _at |
          1  |    56.94289    3.084342     18.46    0.000      50.89769     62.98809
          2  |    54.48244    1.628404     33.46    0.000      51.29082     57.67405
          3  |     52.1283    1.048295     49.73    0.000      50.07368     54.18292
          4  |    49.87588    1.910412     26.11    0.000      46.13154     53.62022
          5  |    47.72079    3.048734     15.65    0.000      41.74538      53.6962
-------------+----------------------------------------------------------------


Variables that uniquely identify margins: milk120


. di exp(3.935802) "  ,  " exp(4.148295) /* correct CI for 1400 */
51.203198  ,  63.325937

.
. * 2a)
. regress wpc i.rp i.vag_disch

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

. margins rp vag_disch

 Predictive margins                                      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 |
         no  |    67.82991    1.366926     49.62    0.000      65.14872     70.51111
        yes  |    78.06962    4.310854     18.11    0.000      69.61398     86.52525
             |
   vag_disch |
         no  |    68.32688     1.33448     51.20    0.000      65.70933     70.94443
        yes  |    77.39382    5.816838     13.31    0.000      65.98424     88.80341
-------------+----------------------------------------------------------------

. table rp vag_disch /* Stata 16: table rp vag_disch, row col */


-----------------------------+--------------------------------
                             |    Vaginal discharge observed
                             |         no       yes      Total
-----------------------------+--------------------------------
Retained placenta at calving |
  no                         |      1,373        52      1,425
  yes                        |        119        30        149
  Total                      |      1,492        82      1,574
-----------------------------+--------------------------------

. lincom _cons+1.rp+1.vag_disch*82/1574 /* rp=1 */


 ( 1)  1.rp + .0520966*1.vag_disch + _cons = 0


-------------+----------------------------------------------------------------
          wpc | Coefficient   Std. err.       t    P>|t|      [95% con f. interval]
-------------+----------------------------------------------------------------
         (1) |    78.06962    4.310854     18.11    0.000      69.61398     86.52525
-------------+----------------------------------------------------------------

. margins rp vag_disch, asbalanced

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


Expression:  Linear prediction, predict()


At:   rp     (asbalanced)
  vag_disch     (asbalanced)

-------------+----------------------------------------------------------------
             |            Delta-method
             |     Margin    std. err.       t    P>|t|      [95% con f. interval]
-------------+----------------------------------------------------------------
          rp |
         no  |    71.89103    3.089892     23.27    0.000      65.83028     77.95177
        yes  |    82.13073    4.580191     17.93    0.000       73.1468     91.11466
             |
   vag_disch |
         no  |    72.47741    2.327566     31.14    0.000      67.91195     77.04287
        yes  |    81.54435     5.71754     14.26    0.000      70.32954     92.75916
-------------+----------------------------------------------------------------

. lincom _cons+1.rp+1.vag_disch*0.5 /* rp=1 */


 ( 1)  1.rp + .5*1.vag_disch + _cons = 0


-------------+----------------------------------------------------------------
          wpc | Coefficient   Std. err.       t    P>|t|      [95% con f. interval]
-------------+----------------------------------------------------------------
         (1) |    82.13073    4.580191     17.93    0.000       73.1468     91.11466
-------------+----------------------------------------------------------------

. margins rp, over(vag_disch) /* same as: at(vag_disch=(0 1)) */

 Predictive margins                                      Number of obs    = 1,574
 Model VCE: OLS


Expression:  Linear prediction, predict()


Over:  vag_disch


-------------+----------------------------------------------------------------
             |            Delta-method
             |     Margin    std. err.       t    P>|t|      [95% con f. interval]
-------------+----------------------------------------------------------------
vag_disch#rp |
      no#no  |    67.35756    1.381095     48.77    0.000      64.64857     70.06654
     no#yes  |    77.59726    4.386003     17.69    0.000      68.99422      86.2003
     yes#no  |     76.4245    5.922753     12.90    0.000      64.80717     88.04183
    yes#yes  |     86.6642    6.372663     13.60    0.000      74.16438     99.16402
-------------+----------------------------------------------------------------

. marginsplot, noci


Variables that uniquely identify margins: rp vag_disch


. lincom _cons+1.rp+1.vag_disch /* rp=1, vag_disch=1 */


 ( 1)  1.rp + 1.vag_disch + _cons = 0


-------------+----------------------------------------------------------------
          wpc | Coefficient   Std. err.       t    P>|t|      [95% con f. interval]
-------------+----------------------------------------------------------------
         (1) |     86.6642    6.372663     13.60    0.000      74.16438     99.16402
-------------+----------------------------------------------------------------

.
. * 2b)
. regress wpc rp##vag_disch

      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 rp#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 /* this is the interaction plot! */


Variables that uniquely identify margins: rp vag_disch


. marginsplot, noci x(vag_disch) /* controlling variable on x */


Variables that uniquely identify margins: rp vag_disch


. marginsplot, noci x(rp) /* same as default */


Variables that uniquely identify margins: rp vag_disch


. marginsplot, noci nosimplelabels /* more informative labels */


Variables that uniquely identify margins: rp vag_disch


. table rp vag_disch, statistic(mean wpc)


-----------------------------+--------------------------------
                             |    Vaginal discharge observed
                             |        no        yes      Total
-----------------------------+--------------------------------
Retained placenta at calving |
  no                         |  67.66861   68.21154   67.68842
  yes                        |   74.0084      100.9   79.42282
  Total                      |  68.17426   80.17073   68.79924
-----------------------------+--------------------------------

. * Stata 16: table rp vag_disch, contents(mean wpc) row col
. margins rp

 Predictive margins                                      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 |
         no  |    67.69689    1.367034     49.52    0.000      65.01549      70.3783
        yes  |    75.40936    4.495364     16.77    0.000      66.59181     84.22691
-------------+----------------------------------------------------------------

. lincom _cons+1.rp+(1.vag_disch+1.rp#1.vag_disch)*82/1574  /* rp=1 */


 ( 1)  1.rp + .0520966*1.vag_disch + .0520966*1.rp#1.vag_disch + _cons = 0


-------------+----------------------------------------------------------------
          wpc | Coefficient   Std. err.       t    P>|t|      [95% con f. interval]
-------------+----------------------------------------------------------------
         (1) |    75.40936    4.495364     16.77    0.000      66.59181     84.22691
-------------+----------------------------------------------------------------

. margins rp, asbalanced

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


Expression:  Linear prediction, predict()


At:   rp     (asbalanced)
  vag_disch     (asbalanced)

-------------+----------------------------------------------------------------
             |            Delta-method
             |     Margin    std. err.       t    P>|t|      [95% con f. interval]
-------------+----------------------------------------------------------------
          rp |
         no  |    67.94007    3.632691     18.70    0.000      60.81464     75.06551
        yes  |     87.4542    5.253116     16.65    0.000      77.15034     97.75806
-------------+----------------------------------------------------------------

. lincom _cons+1.rp+(1.vag_disch+1.rp#1.vag_disch)*0.5  /* rp=1 */


 ( 1)  1.rp + .5*1.vag_disch + .5*1.rp#1.vag_disch + _cons = 0


-------------+----------------------------------------------------------------
          wpc | Coefficient   Std. err.       t    P>|t|      [95% con f. interval]
-------------+----------------------------------------------------------------
         (1) |     87.4542    5.253116     16.65    0.000      77.15034     97.75806
-------------+----------------------------------------------------------------

.
. * 2c)
. 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
-------------+----------------------------------------------------------------

. margins dyst, at( milk120=(1400(1000)5400))

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


Expression:  Linear prediction, predict()


1._at:   milk120   = 1400
2._at:   milk120   = 2400
3._at:   milk120   = 3400
4._at:   milk120   = 4400
5._at:   milk120   = 5400

-------------+----------------------------------------------------------------
             |            Delta-method
             |     Margin    std. err.       t    P>|t|      [95% con f. interval]
-------------+----------------------------------------------------------------
    _at#dyst |
       1#no  |    72.80639    3.713082     19.61    0.000      65.52313     80.08964
      1#yes  |    76.81056    6.244315     12.30    0.000      64.56226     89.05886
       2#no  |    70.56893    2.071173     34.07    0.000       66.5063     74.63156
      2#yes  |    74.57311    5.535781     13.47    0.000       63.7146     85.43161
       3#no  |    68.33148    1.400653     48.79    0.000      65.58408     71.07888
      3#yes  |    72.33565    5.428403     13.33    0.000      61.68777     82.98354
       4#no  |    66.09402    2.608293     25.34    0.000      60.97782     71.21022
      4#yes  |     70.0982    5.954792     11.77    0.000       58.4178      81.7786
       5#no  |    63.85657    4.337493     14.72    0.000      55.34852     72.36462
      5#yes  |    67.86075    6.972863      9.73    0.000      54.18339      81.5381
-------------+----------------------------------------------------------------

. marginsplot, noci


Variables that uniquely identify margins: milk120 dyst


. lincom _cons+1.dyst+milk120*3400 /* dyst=1, milk120=3400 */


 ( 1)  1.dyst + 3400*milk120 + _cons = 0


-------------+----------------------------------------------------------------
          wpc | Coefficient   Std. err.       t    P>|t|      [95% con f. interval]
-------------+----------------------------------------------------------------
         (1) |    72.33565    5.428403     13.33    0.000      61.68777     82.98354
-------------+----------------------------------------------------------------

. margin dyst, atmeans

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


Expression:  Linear prediction, predict()


At:   0.dyst   = .9401042   (mean)
  1.dyst   = .0598958   (mean)
  milk120   = 3215.096   (mean)

-------------+----------------------------------------------------------------
             |            Delta-method
             |     Margin    std. err.       t    P>|t|      [95% con f. interval]
-------------+----------------------------------------------------------------
        dyst |
         no  |    68.74519    1.360888     50.51    0.000      66.07579     71.41459
        yes  |    72.74937     5.39861     13.48    0.000      62.15993     83.33881
-------------+----------------------------------------------------------------

. lincom _cons+1.dyst+milk120*3215.096


 ( 1)  1.dyst + 3215.096*milk120 + _cons = 0


-------------+----------------------------------------------------------------
          wpc | Coefficient   Std. err.       t    P>|t|      [95% con f. interval]
-------------+----------------------------------------------------------------
         (1) |    72.74937     5.39861     13.48    0.000      62.15993     83.33881
-------------+----------------------------------------------------------------

. * interaction model
. regress wpc dyst##c.milk120

      Source |       SS           df       MS      Number of obs   =     1,536
-------------+----------------------------------   F(3, 1532)      =      3.83
       Model |  30572.8757         3  10190.9586   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
       milk120 |   -.0034465    .0019285     -1.79    0.074     -.0072294     .0003363
               |
dyst#c.milk120 |
          yes  |    .0291424    .0094681      3.08    0.002      .0105706     .0477142
               |
         _cons |    79.83774    6.365297     12.54    0.000      67.35213     92.32336
---------------+----------------------------------------------------------------

. margins dyst, at( milk120=(1400(1000)5400))

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


Expression:  Linear prediction, predict()


1._at:   milk120   = 1400
2._at:   milk120   = 2400
3._at:   milk120   = 3400
4._at:   milk120   = 4400
5._at:   milk120   = 5400

-------------+----------------------------------------------------------------
             |            Delta-method
             |     Margin    std. err.       t    P>|t|      [95% con f. interval]
-------------+----------------------------------------------------------------
    _at#dyst |
       1#no  |     75.0126      3.7716     19.89    0.000      67.61456     82.41064
      1#yes  |    30.32354     16.3366      1.86    0.064      -1.72092     62.36801
       2#no  |    71.56607    2.090724     34.23    0.000      67.46508     75.66705
      2#yes  |    56.01939     8.17389      6.85    0.000      39.98619     72.05258
       3#no  |    68.11954    1.398494     48.71    0.000      65.37637      70.8627
      3#yes  |    81.71523    6.212235     13.15    0.000      69.52985     93.90061
       4#no  |      64.673    2.641768     24.48    0.000      59.49114     69.85487
      4#yes  |    107.4111    13.49899      7.96    0.000      80.93263     133.8895
       5#no  |    61.22647    4.409148     13.89    0.000      52.57787     69.87508
      5#yes  |    133.1069     22.3093      5.97    0.000      89.34692     176.8669
-------------+----------------------------------------------------------------

. marginsplot, noci /* this is the interaction plot! */


Variables that uniquely identify margins: milk120 dyst


. lincom _cons+1.dyst+(c.milk120+1.dyst#c.milk120)*3400 /* dyst=1, milk120=3200 */


 ( 1)  1.dyst + 3400*milk120 + 3400*1.dyst#c.milk120 + _cons = 0


-------------+----------------------------------------------------------------
          wpc | Coefficient   Std. err.       t    P>|t|      [95% con f. interval]
-------------+----------------------------------------------------------------
         (1) |    81.71523    6.212235     13.15    0.000      69.52985     93.90061
-------------+----------------------------------------------------------------

.
. * 2d)
. regress wpc parity milk120

      Source |       SS           df       MS      Number of obs   =     1,536
-------------+----------------------------------   F(2, 1533)      =      3.18
       Model |  16982.1525         2  8491.07626   Prob > F        =    0.0417
    Residual |   4087382.5     1,533  2666.26386   R-squared       =    0.0041
-------------+----------------------------------   Adj R-squared   =    0.0028
       Total |  4104364.66     1,535   2673.8532   Root MSE        =    51.636

-------------+----------------------------------------------------------------
          wpc | Coefficient   Std. err.       t    P>|t|      [95% con f. interval]
-------------+----------------------------------------------------------------
      parity |    2.104084    .9534787      2.21    0.027      .2338231     3.974344
     milk120 |   -.0040343    .0020428     -1.97    0.048     -.0080413    -.0000273
       _cons |    76.19823    6.211661     12.27    0.000      64.01398     88.38248
-------------+----------------------------------------------------------------

. tabulate parity

  Lactation |
     number |      Freq.     Percent        Cum.
------------+-----------------------------------
          1 |        417       26.49       26.49
          2 |        374       23.76       50.25
          3 |        319       20.27       70.52
          4 |        222       14.10       84.63
          5 |        169       10.74       95.36
          6 |         69        4.38       99.75
          7 |          4        0.25      100.00
------------+-----------------------------------
      Total |      1,574      100.00

. margins , at( parity=(1(1)6) milk120=(1400(1000)5400))

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


Expression:  Linear prediction, predict()


1._at:   parity   = 1
  milk120   = 1400
2._at:   parity   = 1
  milk120   = 2400
3._at:   parity   = 1
  milk120   = 3400
4._at:   parity   = 1
  milk120   = 4400
5._at:   parity   = 1
  milk120   = 5400
6._at:   parity   = 2
  milk120   = 1400
7._at:   parity   = 2
  milk120   = 2400
8._at:   parity   = 2
  milk120   = 3400
9._at:   parity   = 2
  milk120   = 4400
10._at:   parity   = 2
  milk120   = 5400
11._at:   parity   = 3
  milk120   = 1400
12._at:   parity   = 3
  milk120   = 2400
13._at:   parity   = 3
  milk120   = 3400
14._at:   parity   = 3
  milk120   = 4400
15._at:   parity   = 3
  milk120   = 5400
16._at:   parity   = 4
  milk120   = 1400
17._at:   parity   = 4
  milk120   = 2400
18._at:   parity   = 4
  milk120   = 3400
19._at:   parity   = 4
  milk120   = 4400
20._at:   parity   = 4
  milk120   = 5400
21._at:   parity   = 5
  milk120   = 1400
22._at:   parity   = 5
  milk120   = 2400
23._at:   parity   = 5
  milk120   = 3400
24._at:   parity   = 5
  milk120   = 4400
25._at:   parity   = 5
  milk120   = 5400
26._at:   parity   = 6
  milk120   = 1400
27._at:   parity   = 6
  milk120   = 2400
28._at:   parity   = 6
  milk120   = 3400
29._at:   parity   = 6
  milk120   = 4400
30._at:   parity   = 6
  milk120   = 5400

-------------+----------------------------------------------------------------
             |            Delta-method
             |     Margin    std. err.       t    P>|t|      [95% con f. interval]
-------------+----------------------------------------------------------------
         _at |
          1  |    72.65429    3.678887     19.75    0.000      65.43811     79.87047
          2  |    68.61999    2.267725     30.26    0.000      64.17182     73.06816
          3  |    64.58569    2.257717     28.61    0.000      60.15715     69.01423
          4  |    60.55138    3.660373     16.54    0.000      53.37152     67.73125
          5  |    56.51708    5.481398     10.31    0.000      45.76525     67.26891
          6  |    74.75837    3.740069     19.99    0.000      67.42218     82.09457
          7  |    70.72407    2.026805     34.89    0.000      66.74847     74.69968
          8  |    66.68977      1.6044     41.57    0.000      63.54272     69.83682
          9  |    62.65547    3.063757     20.45    0.000      56.64587     68.66507
         10  |    58.62117    4.954337     11.83    0.000      48.90317     68.33916
         11  |    76.86246    4.032402     19.06    0.000      68.95285     84.77207
         12  |    72.82816    2.211683     32.93    0.000      68.48991      77.1664
         13  |    68.79385    1.367172     50.32    0.000      66.11213     71.47558
         14  |    64.75955    2.681999     24.15    0.000      59.49878     70.02033
         15  |    60.72525    4.567649     13.29    0.000      51.76575     69.68475
         16  |    78.96654    4.511171     17.50    0.000      70.11782     87.81526
         17  |    74.93224    2.737406     27.37    0.000      69.56278      80.3017
         18  |    70.89794    1.726981     41.05    0.000      67.51044     74.28543
         19  |    66.86363    2.611106     25.61    0.000      61.74192     71.98535
         20  |    62.82933    4.358857     14.41    0.000      54.27938     71.37929
         21  |    81.07062    5.124382     15.82    0.000      71.01908     91.12217
         22  |    77.03632    3.451591     22.32    0.000      70.26598     83.80666
         23  |    73.00202    2.431873     30.02    0.000      68.23187     77.77217
         24  |    68.96772    2.874174     24.00    0.000      63.32999     74.60545
         25  |    64.93342    4.353629     14.91    0.000      56.39372     73.47312
         26  |    83.17471    5.829766     14.27    0.000      71.73955     94.60987
         27  |    79.14041    4.260495     18.58    0.000      70.78339     87.49742
         28  |     75.1061    3.265547     23.00    0.000      68.70069     81.51152
         29  |     71.0718    3.394425     20.94    0.000      64.41359     77.73001
         30  |     67.0375    4.552668     14.72    0.000      58.10738     75.96762
-------------+----------------------------------------------------------------

. marginsplot, noci


Variables that uniquely identify margins: parity milk120


. margins , at( milk120=(1400(1000)5400) parity=(1(1)6))

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


Expression:  Linear prediction, predict()


1._at:   parity   = 1
  milk120   = 1400
2._at:   parity   = 1
  milk120   = 2400
3._at:   parity   = 1
  milk120   = 3400
4._at:   parity   = 1
  milk120   = 4400
5._at:   parity   = 1
  milk120   = 5400
6._at:   parity   = 2
  milk120   = 1400
7._at:   parity   = 2
  milk120   = 2400
8._at:   parity   = 2
  milk120   = 3400
9._at:   parity   = 2
  milk120   = 4400
10._at:   parity   = 2
  milk120   = 5400
11._at:   parity   = 3
  milk120   = 1400
12._at:   parity   = 3
  milk120   = 2400
13._at:   parity   = 3
  milk120   = 3400
14._at:   parity   = 3
  milk120   = 4400
15._at:   parity   = 3
  milk120   = 5400
16._at:   parity   = 4
  milk120   = 1400
17._at:   parity   = 4
  milk120   = 2400
18._at:   parity   = 4
  milk120   = 3400
19._at:   parity   = 4
  milk120   = 4400
20._at:   parity   = 4
  milk120   = 5400
21._at:   parity   = 5
  milk120   = 1400
22._at:   parity   = 5
  milk120   = 2400
23._at:   parity   = 5
  milk120   = 3400
24._at:   parity   = 5
  milk120   = 4400
25._at:   parity   = 5
  milk120   = 5400
26._at:   parity   = 6
  milk120   = 1400
27._at:   parity   = 6
  milk120   = 2400
28._at:   parity   = 6
  milk120   = 3400
29._at:   parity   = 6
  milk120   = 4400
30._at:   parity   = 6
  milk120   = 5400

-------------+----------------------------------------------------------------
             |            Delta-method
             |     Margin    std. err.       t    P>|t|      [95% con f. interval]
-------------+----------------------------------------------------------------
         _at |
          1  |    72.65429    3.678887     19.75    0.000      65.43811     79.87047
          2  |    68.61999    2.267725     30.26    0.000      64.17182     73.06816
          3  |    64.58569    2.257717     28.61    0.000      60.15715     69.01423
          4  |    60.55138    3.660373     16.54    0.000      53.37152     67.73125
          5  |    56.51708    5.481398     10.31    0.000      45.76525     67.26891
          6  |    74.75837    3.740069     19.99    0.000      67.42218     82.09457
          7  |    70.72407    2.026805     34.89    0.000      66.74847     74.69968
          8  |    66.68977      1.6044     41.57    0.000      63.54272     69.83682
          9  |    62.65547    3.063757     20.45    0.000      56.64587     68.66507
         10  |    58.62117    4.954337     11.83    0.000      48.90317     68.33916
         11  |    76.86246    4.032402     19.06    0.000      68.95285     84.77207
         12  |    72.82816    2.211683     32.93    0.000      68.48991      77.1664
         13  |    68.79385    1.367172     50.32    0.000      66.11213     71.47558
         14  |    64.75955    2.681999     24.15    0.000      59.49878     70.02033
         15  |    60.72525    4.567649     13.29    0.000      51.76575     69.68475
         16  |    78.96654    4.511171     17.50    0.000      70.11782     87.81526
         17  |    74.93224    2.737406     27.37    0.000      69.56278      80.3017
         18  |    70.89794    1.726981     41.05    0.000      67.51044     74.28543
         19  |    66.86363    2.611106     25.61    0.000      61.74192     71.98535
         20  |    62.82933    4.358857     14.41    0.000      54.27938     71.37929
         21  |    81.07062    5.124382     15.82    0.000      71.01908     91.12217
         22  |    77.03632    3.451591     22.32    0.000      70.26598     83.80666
         23  |    73.00202    2.431873     30.02    0.000      68.23187     77.77217
         24  |    68.96772    2.874174     24.00    0.000      63.32999     74.60545
         25  |    64.93342    4.353629     14.91    0.000      56.39372     73.47312
         26  |    83.17471    5.829766     14.27    0.000      71.73955     94.60987
         27  |    79.14041    4.260495     18.58    0.000      70.78339     87.49742
         28  |     75.1061    3.265547     23.00    0.000      68.70069     81.51152
         29  |     71.0718    3.394425     20.94    0.000      64.41359     77.73001
         30  |     67.0375    4.552668     14.72    0.000      58.10738     75.96762
-------------+----------------------------------------------------------------

. marginsplot, noci /* changing roles in plot */


Variables that uniquely identify margins: milk120 parity


. margins , at( milk120=(1400(1000)5400) (median)parity)

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


Expression:  Linear prediction, predict()


1._at:   parity   = 2   (median)
  milk120   = 1400
2._at:   parity   = 2   (median)
  milk120   = 2400
3._at:   parity   = 2   (median)
  milk120   = 3400
4._at:   parity   = 2   (median)
  milk120   = 4400
5._at:   parity   = 2   (median)
  milk120   = 5400

-------------+----------------------------------------------------------------
             |            Delta-method
             |     Margin    std. err.       t    P>|t|      [95% con f. interval]
-------------+----------------------------------------------------------------
         _at |
          1  |    74.75837    3.740069     19.99    0.000      67.42218     82.09457
          2  |    70.72407    2.026805     34.89    0.000      66.74847     74.69968
          3  |    66.68977      1.6044     41.57    0.000      63.54272     69.83682
          4  |    62.65547    3.063757     20.45    0.000      56.64587     68.66507
          5  |    58.62117    4.954337     11.83    0.000      48.90317     68.33916
-------------+----------------------------------------------------------------

. marginsplot


Variables that uniquely identify margins: milk120


. lincom _cons+milk120*1400+parity*2 /* milk120=1400, parity=2 */


 ( 1)  2*parity + 1400*milk120 + _cons = 0


-------------+----------------------------------------------------------------
          wpc | Coefficient   Std. err.       t    P>|t|      [95% con f. interval]
-------------+----------------------------------------------------------------
         (1) |    74.75837    3.740069     19.99    0.000      67.42218     82.09457
-------------+----------------------------------------------------------------

. margins, atmeans

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


Expression:  Linear prediction, predict()


At:   parity   = 2.736328   (mean)
  milk120   = 3215.096   (mean)

-------------+----------------------------------------------------------------
             |            Delta-method
             |     Margin    std. err.       t    P>|t|      [95% con f. interval]
-------------+----------------------------------------------------------------
       _cons |    68.98503    1.317516     52.36    0.000       66.4007     71.56935
-------------+----------------------------------------------------------------

. lincom _cons+milk120*3215.096+parity*2.736328 /* both at means */


 ( 1)  2.736328*parity + 3215.096*milk120 + _cons = 0


-------------+----------------------------------------------------------------
          wpc | Coefficient   Std. err.       t    P>|t|      [95% con f. interval]
-------------+----------------------------------------------------------------
         (1) |    68.98502    1.317516     52.36    0.000       66.4007     71.56935
-------------+----------------------------------------------------------------

. * interaction model
. regress wpc c.parity##c.milk120

      Source |       SS           df       MS      Number of obs   =     1,536
-------------+----------------------------------   F(3, 1532)      =      2.26
       Model |  18084.1915         3  6028.06382   Prob > F        =    0.0797
    Residual |  4086280.46     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.88999    4.437817      1.10    0.271     -3.814849     13.59483
           milk120 |   -.0016664    .0042125     -0.40    0.692     -.0099293     .0065965
                   |
c.parity#c.milk120 |   -.0008764    .0013635     -0.64    0.520      -.003551     .0017981
                   |
             _cons |    69.02201    12.77659      5.40    0.000      43.96056     94.08347
-------------------+----------------------------------------------------------------

. margins , at( parity=(1(1)6) milk120=(1400(1000)5400))

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


Expression:  Linear prediction, predict()


1._at:   parity   = 1
  milk120   = 1400
2._at:   parity   = 1
  milk120   = 2400
3._at:   parity   = 1
  milk120   = 3400
4._at:   parity   = 1
  milk120   = 4400
5._at:   parity   = 1
  milk120   = 5400
6._at:   parity   = 2
  milk120   = 1400
7._at:   parity   = 2
  milk120   = 2400
8._at:   parity   = 2
  milk120   = 3400
9._at:   parity   = 2
  milk120   = 4400
10._at:   parity   = 2
  milk120   = 5400
11._at:   parity   = 3
  milk120   = 1400
12._at:   parity   = 3
  milk120   = 2400
13._at:   parity   = 3
  milk120   = 3400
14._at:   parity   = 3
  milk120   = 4400
15._at:   parity   = 3
  milk120   = 5400
16._at:   parity   = 4
  milk120   = 1400
17._at:   parity   = 4
  milk120   = 2400
18._at:   parity   = 4
  milk120   = 3400
19._at:   parity   = 4
  milk120   = 4400
20._at:   parity   = 4
  milk120   = 5400
21._at:   parity   = 5
  milk120   = 1400
22._at:   parity   = 5
  milk120   = 2400
23._at:   parity   = 5
  milk120   = 3400
24._at:   parity   = 5
  milk120   = 4400
25._at:   parity   = 5
  milk120   = 5400
26._at:   parity   = 6
  milk120   = 1400
27._at:   parity   = 6
  milk120   = 2400
28._at:   parity   = 6
  milk120   = 3400
29._at:   parity   = 6
  milk120   = 4400
30._at:   parity   = 6
  milk120   = 5400

-------------+----------------------------------------------------------------
             |            Delta-method
             |     Margin    std. err.       t    P>|t|      [95% con f. interval]
-------------+----------------------------------------------------------------
         _at |
          1  |    70.35204    5.134977     13.70    0.000      60.27971     80.42436
          2  |     67.8092    2.595304     26.13    0.000      62.71848     72.89993
          3  |    65.26637     2.49412     26.17    0.000      60.37412     70.15862
          4  |    62.72353    4.982277     12.59    0.000      52.95073     72.49634
          5  |     60.1807     7.90841      7.61    0.000      44.66824     75.69315
          6  |    74.01502    3.915469     18.90    0.000      66.33477     81.69526
          7  |    70.59575       2.037     34.66    0.000      66.60014     74.59135
          8  |    67.17647    1.774378     37.86    0.000      63.69601     70.65694
          9  |     63.7572    3.511131     18.16    0.000      56.87007     70.64434
         10  |    60.33793    5.629231     10.72    0.000      49.29612     71.37975
         11  |    77.67799    4.228031     18.37    0.000      69.38465     85.97134
         12  |    73.38229    2.374154     30.91    0.000      68.72535     78.03922
         13  |    69.08658    1.441275     47.93    0.000       66.2595     71.91366
         14  |    64.79088    2.682955     24.15    0.000      59.52822     70.05353
         15  |    60.49517    4.582525     13.20    0.000      51.50648     69.48385
         16  |    81.34097    5.831298     13.95    0.000       69.9028     92.77914
         17  |    76.16883    3.346239     22.76    0.000      69.60514     82.73253
         18  |    70.99669     1.73413     40.94    0.000      67.59517     74.39821
         19  |    65.82455    3.071435     21.43    0.000      59.79989     71.84921
         20  |     60.6524    5.520583     10.99    0.000      49.82371      71.4811
         21  |    85.00395    7.982122     10.65    0.000      69.34691      100.661
         22  |    78.95537    4.564154     17.30    0.000      70.00273     87.90802
         23  |     72.9068    2.436846     29.92    0.000      68.12689      77.6867
         24  |    66.85822    4.362847     15.32    0.000      58.30043       75.416
         25  |    60.80964    7.753718      7.84    0.000      45.60062     76.01866
         26  |    88.66693    10.34442      8.57    0.000      68.37622     108.9576
         27  |    81.74192    5.877002     13.91    0.000       70.2141     93.26974
         28  |     74.8169    3.297015     22.69    0.000      68.34976     81.28404
         29  |    67.89189    6.000034     11.32    0.000      56.12274     79.66104
         30  |    60.96688    10.48472      5.81    0.000      40.40094     81.53281
-------------+----------------------------------------------------------------

. marginsplot, noci


Variables that uniquely identify margins: parity milk120


. lincom _cons+milk120*1400+parity*2+milk120#c.parity*1400*2 /* milk120=1400, parity=2 */


 ( 1)  2*parity + 1400*milk120 + 2800*c.parity#c.milk120 + _cons = 0


-------------+----------------------------------------------------------------
          wpc | Coefficient   Std. err.       t    P>|t|      [95% con f. interval]
-------------+----------------------------------------------------------------
         (1) |    74.01502    3.915469     18.90    0.000      66.33477     81.69526
-------------+----------------------------------------------------------------

.
. * VER Ex. 14.16: estimates on sqrt-scale and backtransformed
. gen calv_mth=month(calv_dt)

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

. gen rootwpc=sqrt(wpc)

. regress rootwpc aut_calv c.herd_size##c.herd_size parity twin dyst rp##vag_disch

      Source |       SS           df       MS      Number of obs   =     1,574
-------------+----------------------------------   F(9, 1564)      =     16.19
       Model |  1133.93223         9   125.99247   Prob > F        =    0.0000
    Residual |  12172.3346     1,564  7.78282264   R-squared       =    0.0852
-------------+----------------------------------   Adj R-squared   =    0.0800
       Total |  13306.2668     1,573  8.45916519   Root MSE        =    2.7898

------------------------+----------------------------------------------------------------
                 rootwpc | Coefficient   Std. err.       t    P>|t|      [95% con f. interval]
------------------------+----------------------------------------------------------------
               aut_calv |   -.5136529    .1419214     -3.62    0.000     -.7920293    -.2352766
              herd_size |   -.0232678    .0084168     -2.76    0.006     -.0397771    -.0067585
                        |
c.herd_size#c.herd_size |    .0000709    .0000174      4.07    0.000      .0000367      .000105
                        |
                 parity |     .058304    .0480002      1.21    0.225     -.0358476     .1524556
                   twin |    1.385453    .5505819      2.52    0.012      .3054964     2.465409
                   dyst |    .5422828    .3054896      1.78    0.076     -.0569296     1.141495
                        |
                     rp |
                   yes  |    .3894596    .2691054      1.45    0.148     -.1383858      .917305
                        |
              vag_disch |
                   yes  |   -.0132839    .4004945     -0.03    0.974     -.7988467     .7722788
                        |
           rp#vag_disch |
               yes#yes  |    1.491019    .6999486      2.13    0.033      .1180826     2.863956
                        |
                  _cons |    8.834599    .9921073      8.90    0.000      6.888598      10.7806
------------------------+----------------------------------------------------------------

. *1) correct estimate and CI for rp=1, vag_disch=0, dyst=1
. lincom _cons+aut_calv*0+herd_size*234+c.herd_size#c.herd_size*54756+parity*1+twin*0+dyst*1+1.rp+0.vag_disch+1.rp#0.vag_disch


 ( 1)  234*herd_size + 54756*c.herd_size#c.herd_size + parity + dyst + 1.rp + 0b.vag_disch + 1o.rp#0b.vag_disch + _cons = 0


-------------+----------------------------------------------------------------
      rootwpc | Coefficient   Std. err.       t    P>|t|      [95% con f. interval]
-------------+----------------------------------------------------------------
         (1) |    8.259744    .3763213     21.95    0.000      7.521596     8.997891
-------------+----------------------------------------------------------------

. *1) correct estimate and SE using margins command
. margins , over(rp vag_disch dyst) at(parity=1 aut_calv=0 herd_size=234 twin=0)

 Predictive margins                                      Number of obs    = 1,574
 Model VCE: OLS


Expression:  Linear prediction, predict()


Over:  rp vag_disch dyst


At: 0.rp#0.vag_disch#0.dyst
      aut_calv   = 0
      herd_size   = 234
      parity   = 1
      twin   = 0
    0.rp#0.vag_disch#1.dyst
      aut_calv   = 0
      herd_size   = 234
      parity   = 1
      twin   = 0
    0.rp#1.vag_disch#0.dyst
      aut_calv   = 0
      herd_size   = 234
      parity   = 1
      twin   = 0
    0.rp#1.vag_disch#1.dyst
      aut_calv   = 0
      herd_size   = 234
      parity   = 1
      twin   = 0
    1.rp#0.vag_disch#0.dyst
      aut_calv   = 0
      herd_size   = 234
      parity   = 1
      twin   = 0
    1.rp#0.vag_disch#1.dyst
      aut_calv   = 0
      herd_size   = 234
      parity   = 1
      twin   = 0
    1.rp#1.vag_disch#0.dyst
      aut_calv   = 0
      herd_size   = 234
      parity   = 1
      twin   = 0
    1.rp#1.vag_disch#1.dyst
      aut_calv   = 0
      herd_size   = 234
      parity   = 1
      twin   = 0

------------------+----------------------------------------------------------------
                  |            Delta-method
                  |     Margin    std. err.       t    P>|t|      [95% con f. interval]
------------------+----------------------------------------------------------------
rp#vag_disch#dyst |
        no#no#no  |    7.328001    .1498901     48.89    0.000      7.033995     7.622008
       no#no#yes  |    7.870284    .3094335     25.43    0.000      7.263336     8.477233
       no#yes#no  |    7.314717     .403061     18.15    0.000      6.524121     8.105314
      no#yes#yes  |       7.857    .4546238     17.28    0.000      6.965264     8.748737
       yes#no#no  |    7.717461    .2812858     27.44    0.000      7.165724     8.269198
      yes#no#yes  |    8.259744    .3763213     21.95    0.000      7.521596     8.997891
      yes#yes#no  |    9.195196     .529305     17.37    0.000      8.156974     10.23342
     yes#yes#yes  |    9.737479     .591966     16.45    0.000      8.576348     10.89861
------------------+----------------------------------------------------------------

. *1) correct backtransformed estimate, but inappropriate CI
. margins , over(rp vag_disch dyst) at(parity=1 aut_calv=0 herd_size=234 twin=0) expr(predict(xb)^2)

 Predictive margins                                      Number of obs    = 1,574
 Model VCE: OLS


Expression:  predict(xb)^2


Over:  rp vag_disch dyst


At: 0.rp#0.vag_disch#0.dyst
      aut_calv   = 0
      herd_size   = 234
      parity   = 1
      twin   = 0
    0.rp#0.vag_disch#1.dyst
      aut_calv   = 0
      herd_size   = 234
      parity   = 1
      twin   = 0
    0.rp#1.vag_disch#0.dyst
      aut_calv   = 0
      herd_size   = 234
      parity   = 1
      twin   = 0
    0.rp#1.vag_disch#1.dyst
      aut_calv   = 0
      herd_size   = 234
      parity   = 1
      twin   = 0
    1.rp#0.vag_disch#0.dyst
      aut_calv   = 0
      herd_size   = 234
      parity   = 1
      twin   = 0
    1.rp#0.vag_disch#1.dyst
      aut_calv   = 0
      herd_size   = 234
      parity   = 1
      twin   = 0
    1.rp#1.vag_disch#0.dyst
      aut_calv   = 0
      herd_size   = 234
      parity   = 1
      twin   = 0
    1.rp#1.vag_disch#1.dyst
      aut_calv   = 0
      herd_size   = 234
      parity   = 1
      twin   = 0

------------------+----------------------------------------------------------------
                  |            Delta-method
                  |     Margin    std. err.       z    P>|z|      [95% con f. interval]
------------------+----------------------------------------------------------------
rp#vag_disch#dyst |
        no#no#no  |     53.6996    2.196789     24.44    0.000      49.39398     58.00523
       no#no#yes  |    61.94137     4.87066     12.72    0.000      52.39506     71.48769
       no#yes#no  |    53.50509    5.896555      9.07    0.000      41.94806     65.06213
      no#yes#yes  |    61.73245    7.143958      8.64    0.000      47.73055     75.73435
       yes#no#no  |     59.5592    4.341625     13.72    0.000      51.04978     68.06863
      yes#no#yes  |    68.22337    6.216636     10.97    0.000      56.03899     80.40775
      yes#yes#no  |    84.55163    9.734126      8.69    0.000       65.4731     103.6302
     yes#yes#yes  |     94.8185    11.52851      8.22    0.000      72.22303      117.414
------------------+----------------------------------------------------------------

. table (rp) (vag_disch dyst), statistic(count wpc) nototals


-----------------------------+-----------------------------------------------
                             |            Vaginal discharge observed
                             |            no                     yes
                             |   Dystocia at calving     Dystocia at calving
                             |          no        yes          no         yes
-----------------------------+-----------------------------------------------
Retained placenta at calving |
  no                         |       1,307         66          40          12
  yes                        |         105         14          28           2
-----------------------------+-----------------------------------------------

. * Stata 16: table rp dyst vag_disch, contents(count wpc)
. table parity


-----------------+-----------
                 |  Frequency
-----------------+-----------
Lactation number |
  1              |        417
  2              |        374
  3              |        319
  4              |        222
  5              |        169
  6              |         69
  7              |          4
  Total          |      1,574
-----------------+-----------

. table twin


-----------+-----------
           |  Frequency
-----------+-----------
Twins born |
  no       |      1,547
  yes      |         27
  Total    |      1,574
-----------+-----------

. table herd_size


----------+-----------
          |  Frequency
----------+-----------
Herd size |
  125     |        129
  185     |        168
  201     |        203
  235     |        227
  263     |        253
  294     |        272
  333     |        322
  Total   |      1,574
----------+-----------

. table aut_calv


---------+-----------
         |  Frequency
---------+-----------
aut_calv |
  0      |        831
  1      |        743
  Total  |      1,574
---------+-----------

. *3) correct estimate and SE using margins command
. margins , at(parity=1 aut_calv=0 herd_size=(125 185 201 235 263 294 333) twin=0 rp=0 dyst=0 vag_disch=0)

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


Expression:  Linear prediction, predict()


1._at:   aut_calv   = 0
  herd_size   = 125
  parity   = 1
  twin   = 0
  dyst   = 0
  rp   = 0
  vag_disch   = 0
2._at:   aut_calv   = 0
  herd_size   = 185
  parity   = 1
  twin   = 0
  dyst   = 0
  rp   = 0
  vag_disch   = 0
3._at:   aut_calv   = 0
  herd_size   = 201
  parity   = 1
  twin   = 0
  dyst   = 0
  rp   = 0
  vag_disch   = 0
4._at:   aut_calv   = 0
  herd_size   = 235
  parity   = 1
  twin   = 0
  dyst   = 0
  rp   = 0
  vag_disch   = 0
5._at:   aut_calv   = 0
  herd_size   = 263
  parity   = 1
  twin   = 0
  dyst   = 0
  rp   = 0
  vag_disch   = 0
6._at:   aut_calv   = 0
  herd_size   = 294
  parity   = 1
  twin   = 0
  dyst   = 0
  rp   = 0
  vag_disch   = 0
7._at:   aut_calv   = 0
  herd_size   = 333
  parity   = 1
  twin   = 0
  dyst   = 0
  rp   = 0
  vag_disch   = 0

-------------+----------------------------------------------------------------
             |            Delta-method
             |     Margin    std. err.       t    P>|t|      [95% con f. interval]
-------------+----------------------------------------------------------------
         _at |
          1  |    7.091544    .2559118     27.71    0.000      6.589578      7.59351
          2  |    7.013389    .1550613     45.23    0.000      6.709239     7.317539
          3  |    7.078708    .1514358     46.74    0.000      6.781669     7.375746
          4  |    7.337965    .1497933     48.99    0.000      7.044148     7.631782
          5  |    7.674476    .1445034     53.11    0.000      7.391035     7.957917
          6  |    8.176637    .1404317     58.23    0.000      7.901183     8.452091
          7  |    9.001823    .1805099     49.87    0.000      8.647756      9.35589
-------------+----------------------------------------------------------------

. *3) correct backtransformed estimate, but inappropriate CI
. margins , at(parity=1 aut_calv=0 herd_size=(125 185 201 235 263 294 333) twin=0 rp=0 dyst=0 vag_disch=0) expr(predict(xb)^2) plot

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


Expression:  predict(xb)^2


1._at:   aut_calv   = 0
  herd_size   = 125
  parity   = 1
  twin   = 0
  dyst   = 0
  rp   = 0
  vag_disch   = 0
2._at:   aut_calv   = 0
  herd_size   = 185
  parity   = 1
  twin   = 0
  dyst   = 0
  rp   = 0
  vag_disch   = 0
3._at:   aut_calv   = 0
  herd_size   = 201
  parity   = 1
  twin   = 0
  dyst   = 0
  rp   = 0
  vag_disch   = 0
4._at:   aut_calv   = 0
  herd_size   = 235
  parity   = 1
  twin   = 0
  dyst   = 0
  rp   = 0
  vag_disch   = 0
5._at:   aut_calv   = 0
  herd_size   = 263
  parity   = 1
  twin   = 0
  dyst   = 0
  rp   = 0
  vag_disch   = 0
6._at:   aut_calv   = 0
  herd_size   = 294
  parity   = 1
  twin   = 0
  dyst   = 0
  rp   = 0
  vag_disch   = 0
7._at:   aut_calv   = 0
  herd_size   = 333
  parity   = 1
  twin   = 0
  dyst   = 0
  rp   = 0
  vag_disch   = 0

-------------+----------------------------------------------------------------
             |            Delta-method
             |     Margin    std. err.       z    P>|z|      [95% con f. interval]
-------------+----------------------------------------------------------------
         _at |
          1  |       50.29    3.629619     13.86    0.000      43.17607     57.40392
          2  |    49.18762     2.17501     22.61    0.000      44.92468     53.45056
          3  |     50.1081    2.143939     23.37    0.000      45.90606     54.31015
          4  |    53.84573    2.198356     24.49    0.000      49.53703     58.15443
          5  |    58.89758    2.217976     26.55    0.000      54.55043     63.24473
          6  |    66.85739    2.296517     29.11    0.000       62.3563     71.35848
          7  |    81.03282    3.249837     24.93    0.000      74.66326     87.40238
-------------+----------------------------------------------------------------


Variables that uniquely identify margins: herd_size


.
.
end of do-file

. exit, clear