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

. * do-file for lecture 12a of VHM 812/802, Winter 2022
. version 17 /* works also with versions 14-16 */

. set more off

. cd "r:\"
r:\

.
. * unadjusted and mixed model results (for reference)
. use simcont_clustherd.dta, clear

. regress milk X

      Source |       SS           df       MS      Number of obs   =    11,626
-------------+----------------------------------   F(1, 11624)     =    317.72
       Model |  36598.5078         1  36598.5078   Prob > F        =    0.0000
    Residual |     1338999    11,624  115.192618   R-squared       =    0.0266
-------------+----------------------------------   Adj R-squared   =    0.0265
       Total |  1375597.51    11,625  118.330968   Root MSE        =    10.733

-------------+----------------------------------------------------------------
         milk | Coefficient   Std. err.       t    P>|t|      [95% con f. interval]
-------------+----------------------------------------------------------------
           X |     3.55661     .199534     17.82    0.000       3.16549      3.94773
       _cons |     30.0215    .1457715    205.95    0.000      29.73576     30.30723
-------------+----------------------------------------------------------------

. mixed milk X || herd:, reml

Performing EM optimization ...

Performing gradient-based optimization:
Iteration 0:   log restricted-likelihood = -40902.479
Iteration 1:   log restricted-likelihood = -40902.479  (backed up)

Computing standard errors ...

Mixed-effects REML regression                   Number of obs     =     11,626
Group variable: herd                            Number of groups  =        100
                                                Obs per group:
                                                              min =         20
                                                              avg =      116.3
                                                              max =        311
                                                Wald chi2(1)      =       6.44
Log restricted-likelihood = -40902.479          Prob > chi2       =     0.0112

-------------+----------------------------------------------------------------
         milk | Coefficient   Std. err.       z    P>|z|      [95% con f. interval]
-------------+----------------------------------------------------------------
           X |    3.796004    1.495942      2.54    0.011      .8640104     6.727997
       _cons |    31.13696    1.058717     29.41    0.000      29.06191     33.21201
-------------+----------------------------------------------------------------

-----------------------------+------------------------------------------------
  Random-effects parameters  |   Estimate   Std. err.     [95% conf. interval]
-----------------------------+------------------------------------------------
herd: Identity               |
                  var(_cons) |   54.91494   7.998609       41.2771    73.05868
-----------------------------+------------------------------------------------
               var(Residual) |   64.20087   .8457062      62.56453    65.88001
-----------------------------+------------------------------------------------
LR test vs. linear model: chibar2(01) = 6374.40       Prob >= chibar2 =  0.0000

. collapse milk X, by(herd)

. regress milk X

      Source |       SS           df       MS      Number of obs   =       100
-------------+----------------------------------   F(1, 98)        =      6.37
       Model |   356.97798         1   356.97798   Prob > F        =    0.0132
    Residual |  5493.56229        98   56.056758   R-squared       =    0.0610
-------------+----------------------------------   Adj R-squared   =    0.0514
       Total |  5850.54027        99  59.0963664   Root MSE        =    7.4871

-------------+----------------------------------------------------------------
         milk | Coefficient   Std. err.       t    P>|t|      [95% con f. interval]
-------------+----------------------------------------------------------------
           X |    3.778772    1.497421      2.52    0.013      .8071885     6.750356
       _cons |    31.16586    1.058837     29.43    0.000      29.06463     33.26708
-------------+----------------------------------------------------------------

. use simcont_clustcow.dta, clear

. regress milk X

      Source |       SS           df       MS      Number of obs   =    11,626
-------------+----------------------------------   F(1, 11624)     =    624.90
       Model |  72138.7619         1  72138.7619   Prob > F        =    0.0000
    Residual |  1341880.62    11,624  115.440522   R-squared       =    0.0510
-------------+----------------------------------   Adj R-squared   =    0.0509
       Total |  1414019.39    11,625  121.636076   Root MSE        =    10.744

-------------+----------------------------------------------------------------
         milk | Coefficient   Std. err.       t    P>|t|      [95% con f. interval]
-------------+----------------------------------------------------------------
           X |    4.982006    .1992962     25.00    0.000      4.591352      5.37266
       _cons |    29.25664    .1412627    207.11    0.000      28.97974     29.53354
-------------+----------------------------------------------------------------

. mixed milk X || herd:, reml

Performing EM optimization ...

Performing gradient-based optimization:
Iteration 0:   log restricted-likelihood = -40947.175
Iteration 1:   log restricted-likelihood = -40947.175  (backed up)

Computing standard errors ...

Mixed-effects REML regression                   Number of obs     =     11,626
Group variable: herd                            Number of groups  =        100
                                                Obs per group:
                                                              min =         20
                                                              avg =      116.3
                                                              max =        311
                                                Wald chi2(1)      =    1108.56
Log restricted-likelihood = -40947.175          Prob > chi2       =     0.0000

-------------+----------------------------------------------------------------
         milk | Coefficient   Std. err.       z    P>|z|      [95% con f. interval]
-------------+----------------------------------------------------------------
           X |    4.968194    .1492174     33.30    0.000      4.675733     5.260655
       _cons |    30.64647    .7281274     42.09    0.000      29.21936     32.07357
-------------+----------------------------------------------------------------

-----------------------------+------------------------------------------------
  Random-effects parameters  |   Estimate   Std. err.     [95% conf. interval]
-----------------------------+------------------------------------------------
herd: Identity               |
                  var(_cons) |   51.41187   7.459585      38.68643    68.32319
-----------------------------+------------------------------------------------
               var(Residual) |   64.71069   .8524578      63.06129    66.40324
-----------------------------+------------------------------------------------
LR test vs. linear model: chibar2(01) = 6310.00       Prob >= chibar2 =  0.0000

. use simbin_clustherd.dta, clear

. logit Y X

Iteration 0:   log likelihood = -6894.3552
Iteration 1:   log likelihood = -6815.0583
Iteration 2:   log likelihood = -6814.7785
Iteration 3:   log likelihood = -6814.7785

Logistic regression                                     Number of obs    = 11,626
                                                        LR chi2(1)    = 159.15
                                                        Prob > chi2    = 0.0000
Log likelihood = -6814.7785                             Pseudo R2     = 0.0115

-------------+----------------------------------------------------------------
            Y | Coefficient   Std. err.       z    P>|z|      [95% con f. interval]
-------------+----------------------------------------------------------------
           X |    .5287317    .0423191     12.49    0.000      .4457877     .6116757
       _cons |   -1.241768    .0325699    -38.13    0.000     -1.305604    -1.177932
-------------+----------------------------------------------------------------

. melogit Y X || herd:

Fitting fixed-effects model:

Iteration 0:   log likelihood = -6828.9777
Iteration 1:   log likelihood = -6814.7876
Iteration 2:   log likelihood = -6814.7785
Iteration 3:   log likelihood = -6814.7785

Refining starting values:

Grid node 0:   log likelihood =  -6065.269

Fitting full model:

Iteration 0:   log likelihood =  -6065.269
Iteration 1:   log likelihood =  -6065.089
Iteration 2:   log likelihood = -6065.0864
Iteration 3:   log likelihood = -6065.0864

Mixed-effects logistic regression               Number of obs     =     11,626
Group variable: herd                            Number of groups  =        100

                                                Obs per group:
                                                              min =         20
                                                              avg =      116.3
                                                              max =        311

Integration method: mvaghermite                 Integration pts.  =          7

                                                Wald chi2(1)      =       9.26
Log likelihood = -6065.0864                     Prob > chi2       =     0.0023
-------------+----------------------------------------------------------------
            Y | Coefficient   Std. err.       z    P>|z|      [95% con f. interval]
-------------+----------------------------------------------------------------
           X |    .6199967    .2037578      3.04    0.002      .2206389     1.019355
       _cons |   -1.305448    .1454551     -8.97    0.000     -1.590534    -1.020361
-------------+----------------------------------------------------------------
 herd         |
   var(_cons)|    .9417563    .1493109                    .6902154     1.284968
-------------+----------------------------------------------------------------
LR test vs. logistic model: chibar2(01) = 1499.38     Prob >= chibar2 =  0.0000

. use simbin_clustcow.dta, clear

. logit Y X

Iteration 0:   log likelihood = -6910.3442
Iteration 1:   log likelihood =   -6811.48
Iteration 2:   log likelihood = -6811.0741
Iteration 3:   log likelihood = -6811.0741

Logistic regression                                     Number of obs    = 11,626
                                                        LR chi2(1)    = 198.54
                                                        Prob > chi2    = 0.0000
Log likelihood = -6811.0741                             Pseudo R2     = 0.0144

-------------+----------------------------------------------------------------
            Y | Coefficient   Std. err.       z    P>|z|      [95% con f. interval]
-------------+----------------------------------------------------------------
           X |    .5863084    .0419748     13.97    0.000      .5040393     .6685775
       _cons |    -1.25032    .0316033    -39.56    0.000     -1.312261    -1.188379
-------------+----------------------------------------------------------------

. melogit Y X || herd:

Fitting fixed-effects model:

Iteration 0:   log likelihood =  -6824.417
Iteration 1:   log likelihood = -6811.0819
Iteration 2:   log likelihood = -6811.0741
Iteration 3:   log likelihood = -6811.0741

Refining starting values:

Grid node 0:   log likelihood = -5999.0535

Fitting full model:

Iteration 0:   log likelihood = -5999.0535
Iteration 1:   log likelihood = -5995.9716
Iteration 2:   log likelihood = -5995.9694
Iteration 3:   log likelihood = -5995.9694

Mixed-effects logistic regression               Number of obs     =     11,626
Group variable: herd                            Number of groups  =        100

                                                Obs per group:
                                                              min =         20
                                                              avg =      116.3
                                                              max =        311

Integration method: mvaghermite                 Integration pts.  =          7

                                                Wald chi2(1)      =     229.28
Log likelihood = -5995.9694                     Prob > chi2       =     0.0000
-------------+----------------------------------------------------------------
            Y | Coefficient   Std. err.       z    P>|z|      [95% con f. interval]
-------------+----------------------------------------------------------------
           X |    .6974798     .046063     15.14    0.000      .6071979     .7877616
       _cons |   -1.361196    .1111563    -12.25    0.000     -1.579059    -1.143334
-------------+----------------------------------------------------------------
 herd         |
   var(_cons)|    1.068314    .1682536                    .7845836      1.45465
-------------+----------------------------------------------------------------
LR test vs. logistic model: chibar2(01) = 1630.21     Prob >= chibar2 =  0.0000

.
. * 2-level somatic cell count data
. use scc40_2level, clear

. mixed t_lnscc h_size c_heifer i.t_season t_dim || herdid:, reml

Performing EM optimization ...

Performing gradient-based optimization:
Iteration 0:   log restricted-likelihood = -3624.9622
Iteration 1:   log restricted-likelihood = -3624.9622

Computing standard errors ...

Mixed-effects REML regression                   Number of obs     =      2,178
Group variable: herdid                          Number of groups  =         40
                                                Obs per group:
                                                              min =         12
                                                              avg =       54.5
                                                              max =        105
                                                Wald chi2(6)      =     244.36
Log restricted-likelihood = -3624.9622          Prob > chi2       =     0.0000

-------------+----------------------------------------------------------------
      t_lnscc | Coefficient   Std. err.       z    P>|z|      [95% con f. interval]
-------------+----------------------------------------------------------------
      h_size |    .0040837    .0037726      1.08    0.279     -.0033105     .0114778
    c_heifer |   -.7367168    .0554447    -13.29    0.000     -.8453863    -.6280472
             |
    t_season |
    apr-jun  |    .1609431    .0906574      1.78    0.076     -.0167422     .3386285
    jul-sep  |    .0015031    .0863774      0.02    0.986     -.1677935     .1707997
    oct-dec  |    .0014582    .0918936      0.02    0.987     -.1786499     .1815663
             |
       t_dim |    .0027731    .0004991      5.56    0.000      .0017949     .0037513
       _cons |    4.641202    .1974215     23.51    0.000      4.254263     5.028141
-------------+----------------------------------------------------------------

-----------------------------+------------------------------------------------
  Random-effects parameters  |   Estimate   Std. err.     [95% conf. interval]
-----------------------------+------------------------------------------------
herdid: Identity             |
                  var(_cons) |   .1491533   .0436191      .0840821    .2645832
-----------------------------+------------------------------------------------
               var(Residual) |   1.557228   .0477206      1.466451    1.653625
-----------------------------+------------------------------------------------
LR test vs. linear model: chibar2(01) = 97.01         Prob >= chibar2 =  0.0000

. mixed t_lnscc h_size c_heifer i.t_season t_dim || herdid:, reml stddev

Performing EM optimization ...

Performing gradient-based optimization:
Iteration 0:   log restricted-likelihood = -3624.9622
Iteration 1:   log restricted-likelihood = -3624.9622

Computing standard errors ...

Mixed-effects REML regression                   Number of obs     =      2,178
Group variable: herdid                          Number of groups  =         40
                                                Obs per group:
                                                              min =         12
                                                              avg =       54.5
                                                              max =        105
                                                Wald chi2(6)      =     244.36
Log restricted-likelihood = -3624.9622          Prob > chi2       =     0.0000

-------------+----------------------------------------------------------------
      t_lnscc | Coefficient   Std. err.       z    P>|z|      [95% con f. interval]
-------------+----------------------------------------------------------------
      h_size |    .0040837    .0037726      1.08    0.279     -.0033105     .0114778
    c_heifer |   -.7367168    .0554447    -13.29    0.000     -.8453863    -.6280472
             |
    t_season |
    apr-jun  |    .1609431    .0906574      1.78    0.076     -.0167422     .3386285
    jul-sep  |    .0015031    .0863774      0.02    0.986     -.1677935     .1707997
    oct-dec  |    .0014582    .0918936      0.02    0.987     -.1786499     .1815663
             |
       t_dim |    .0027731    .0004991      5.56    0.000      .0017949     .0037513
       _cons |    4.641202    .1974215     23.51    0.000      4.254263     5.028141
-------------+----------------------------------------------------------------

-----------------------------+------------------------------------------------
  Random-effects parameters  |   Estimate   Std. err.     [95% conf. interval]
-----------------------------+------------------------------------------------
herdid: Identity             |
                   sd(_cons) |   .3862037   .0564716      .2899691    .5143765
-----------------------------+------------------------------------------------
                sd(Residual) |   1.247889   .0191205      1.210971    1.285933
-----------------------------+------------------------------------------------
LR test vs. linear model: chibar2(01) = 97.01         Prob >= chibar2 =  0.0000

. di 1.96*sqrt(.1491533)
.75695926

. di .1491533/(.1491533+1.557228)  /* ICC and VPC */
.08740913

. testparm i.t_season  /* Wald test for season */


 ( 1)  [t_lnscc]2.t_season = 0


 ( 2)  [t_lnscc]3.t_season = 0


 ( 3)  [t_lnscc]4.t_season = 0


           chi2(  3) =    6.21
         Prob > chi2 =    0.1017

. estimates store full

. mixed t_lnscc h_size c_heifer i.t_season t_dim, reml

Mixed-effects REML regression                   Number of obs     =      2,178
                                                Wald chi2(6)      =     250.46
Log restricted-likelihood = -3673.4664          Prob > chi2       =     0.0000

-------------+----------------------------------------------------------------
      t_lnscc | Coefficient   Std. err.       z    P>|z|      [95% con f. interval]
-------------+----------------------------------------------------------------
      h_size |    .0073914     .001605      4.61    0.000      .0042456     .0105371
    c_heifer |   -.7258892    .0567071    -12.80    0.000     -.8370332    -.6147453
             |
    t_season |
    apr-jun  |    .1058737    .0903724      1.17    0.241     -.0712529     .2830004
    jul-sep  |   -.0454035    .0876933     -0.52    0.605     -.2172792     .1264723
    oct-dec  |   -.0385968     .093818     -0.41    0.681     -.2224768     .1452832
             |
       t_dim |    .0030497    .0005049      6.04    0.000      .0020601     .0040393
       _cons |    4.492766    .1067432     42.09    0.000      4.283553     4.701979
-------------+----------------------------------------------------------------

-----------------------------+------------------------------------------------
  Random-effects parameters  |   Estimate   Std. err.     [95% conf. interval]
-----------------------------+------------------------------------------------
               var(Residual) |   1.679409   .0508912      1.582569    1.782176
-----------------------------+------------------------------------------------

. lrtest full

Likelihood-ratio test
Assumption: . nested within full

LR chi2(1)   = 97.01
Prob > chi2   = 0.0000


Note: The reported degrees of freedom assumes the null hypothesis is not on the boundary of the parameter space. If this is not true, then the reported test is conservative.


Note: LR tests based on REML are valid
only when the fixed-effects specification is identical for
both models.


. mixed t_lnscc || herdid:, reml /* "null" model ~ no predictors */

Performing EM optimization ...

Performing gradient-based optimization:
Iteration 0:   log restricted-likelihood = -3722.5104
Iteration 1:   log restricted-likelihood = -3722.5104

Computing standard errors ...

Mixed-effects REML regression                   Number of obs     =      2,178
Group variable: herdid                          Number of groups  =         40
                                                Obs per group:
                                                              min =         12
                                                              avg =       54.5
                                                              max =        105
                                                Wald chi2(0)      =          .
Log restricted-likelihood = -3722.5104          Prob > chi2       =          .

-------------+----------------------------------------------------------------
      t_lnscc | Coefficient   Std. err.       z    P>|z|      [95% con f. interval]
-------------+----------------------------------------------------------------
       _cons |    4.746742    .0679798     69.83    0.000      4.613505      4.87998
-------------+----------------------------------------------------------------

-----------------------------+------------------------------------------------
  Random-effects parameters  |   Estimate   Std. err.     [95% conf. interval]
-----------------------------+------------------------------------------------
herdid: Identity             |
                  var(_cons) |   .1482566   .0426236      .0843907    .2604553
-----------------------------+------------------------------------------------
               var(Residual) |   1.730413   .0529421      1.629699    1.837352
-----------------------------+------------------------------------------------
LR test vs. linear model: chibar2(01) = 101.05        Prob >= chibar2 =  0.0000

.
. * Reunion Island data, 3-level model for (natural) log cfs
. use reu_cfs, clear

. mixed lncfs i.heifer || herd: || cow:, reml

Performing EM optimization ...

Performing gradient-based optimization:
Iteration 0:   log restricted-likelihood = -1500.2224
Iteration 1:   log restricted-likelihood = -1495.8851
Iteration 2:   log restricted-likelihood = -1495.8804
Iteration 3:   log restricted-likelihood = -1495.8804

Computing standard errors ...

Mixed-effects REML regression                   Number of obs     =      3,027

        Grouping information
        ----------------+--------------------------------------------
                        |     No. of       Observations per group
         Group variable |     groups    Minimum    Average    Maximum
        ----------------+--------------------------------------------
                   herd |         50         13       60.5        226
                    cow |      1,575          1        1.9          5
        ----------------+--------------------------------------------

                                                Wald chi2(1)      =       0.12
Log restricted-likelihood = -1495.8804          Prob > chi2       =     0.7341

-------------+----------------------------------------------------------------
        lncfs | Coefficient   Std. err.       z    P>|z|      [95% con f. interval]
-------------+----------------------------------------------------------------
      heifer |
primiparous  |   -.0058409     .017194     -0.34    0.734     -.0395405     .0278587
       _cons |     4.21903    .0195052    216.30    0.000      4.180801      4.25726
-------------+----------------------------------------------------------------

-----------------------------+------------------------------------------------
  Random-effects parameters  |   Estimate   Std. err.     [95% conf. interval]
-----------------------------+------------------------------------------------
herd: Identity               |
                  var(_cons) |   .0145193   .0036025      .0089279    .0236125
-----------------------------+------------------------------------------------
cow: Identity                |
                  var(_cons) |   .0200572   .0040814      .0134605    .0298868
-----------------------------+------------------------------------------------
               var(Residual) |   .1341146   .0048118      .1250077     .143885
-----------------------------+------------------------------------------------
LR test vs. linear model: chi2(2) = 221.39                Prob > chi2 = 0.0000


Note: LR test is conservative and provided only for reference.


. * NOT: mixed lncfs i.heifer || cow: || herd:, reml
. * ICC: same herd
. di .0145193/(.0145193+.0200572+.1341146)
.08607034

. * ICC: same cow
. di (.0145193+.0200572)/(.0145193+.0200572+.1341146)
.20496932

. * ICC calculation by Stata
. estat icc

Residual intraclass correlation

 -----------------------------+------------------------------------------------
                        Level |        ICC   Std. err.     [95% conf. interval]
 -----------------------------+------------------------------------------------
                         herd |   .0860701   .0196752     .0545356     .1332689
                     cow|herd |   .2049693   .0285332     .1546347     .2665226
 -----------------------------+------------------------------------------------

. * 4-level model with region random effects
. mixed lncfs i.heifer || region: || herd: || cow:, reml

Performing EM optimization ...

Performing gradient-based optimization:
Iteration 0:   log restricted-likelihood =  -1500.237
Iteration 1:   log restricted-likelihood = -1495.8809
Iteration 2:   log restricted-likelihood = -1495.8694
Iteration 3:   log restricted-likelihood = -1495.8694

Computing standard errors ...

Mixed-effects REML regression                   Number of obs     =      3,027

        Grouping information
        ----------------+--------------------------------------------
                        |     No. of       Observations per group
         Group variable |     groups    Minimum    Average    Maximum
        ----------------+--------------------------------------------
                 region |          5        193      605.4        960
                   herd |         50         13       60.5        226
                    cow |      1,575          1        1.9          5
        ----------------+--------------------------------------------

                                                Wald chi2(1)      =       0.12
Log restricted-likelihood = -1495.8694          Prob > chi2       =     0.7308

-------------+----------------------------------------------------------------
        lncfs | Coefficient   Std. err.       z    P>|z|      [95% con f. interval]
-------------+----------------------------------------------------------------
      heifer |
primiparous  |   -.0059168    .0171944     -0.34    0.731     -.0396173     .0277836
       _cons |    4.218036    .0206238    204.52    0.000      4.177614     4.258458
-------------+----------------------------------------------------------------

-----------------------------+------------------------------------------------
  Random-effects parameters  |   Estimate   Std. err.     [95% conf. interval]
-----------------------------+------------------------------------------------
region: Identity             |
                  var(_cons) |   .0001844   .0013431      1.16e-10    292.2396
-----------------------------+------------------------------------------------
herd: Identity               |
                  var(_cons) |   .0144022   .0036567       .008756    .0236892
-----------------------------+------------------------------------------------
cow: Identity                |
                  var(_cons) |   .0200442   .0040818      .0134478    .0298764
-----------------------------+------------------------------------------------
               var(Residual) |   .1341222   .0048124      .1250142    .1438939
-----------------------------+------------------------------------------------
LR test vs. linear model: chi2(3) = 221.41                Prob > chi2 = 0.0000


Note: LR test is conservative and provided only for reference.


. * 4-level model with region fixed effects
. mixed lncfs i.heifer i.region || herd: || cow:, reml

Performing EM optimization ...

Performing gradient-based optimization:
Iteration 0:   log restricted-likelihood = -1505.8762
Iteration 1:   log restricted-likelihood =  -1501.496
Iteration 2:   log restricted-likelihood = -1501.4927
Iteration 3:   log restricted-likelihood = -1501.4927

Computing standard errors ...

Mixed-effects REML regression                   Number of obs     =      3,027

        Grouping information
        ----------------+--------------------------------------------
                        |     No. of       Observations per group
         Group variable |     groups    Minimum    Average    Maximum
        ----------------+--------------------------------------------
                   herd |         50         13       60.5        226
                    cow |      1,575          1        1.9          5
        ----------------+--------------------------------------------

                                                Wald chi2(5)      =       4.10
Log restricted-likelihood = -1501.4927          Prob > chi2       =     0.5347

-------------+----------------------------------------------------------------
        lncfs | Coefficient   Std. err.       z    P>|z|      [95% con f. interval]
-------------+----------------------------------------------------------------
      heifer |
primiparous  |   -.0065108    .0172027     -0.38    0.705     -.0402274     .0272058
             |
      region |
          2  |    .0673308    .0854623      0.79    0.431     -.1001721     .2348338
          3  |    .0578045    .1019399      0.57    0.571      -.141994     .2576029
          4  |    .1382983    .0832458      1.66    0.097     -.0248604      .301457
          5  |    .0780979     .082739      0.94    0.345     -.0840675     .2402634
             |
       _cons |    4.131972     .075647     54.62    0.000      3.983706     4.280237
-------------+----------------------------------------------------------------

-----------------------------+------------------------------------------------
  Random-effects parameters  |   Estimate   Std. err.     [95% conf. interval]
-----------------------------+------------------------------------------------
herd: Identity               |
                  var(_cons) |   .0146174   .0037663      .0088217    .0242209
-----------------------------+------------------------------------------------
cow: Identity                |
                  var(_cons) |   .0199372   .0040785      .0133517    .0297709
-----------------------------+------------------------------------------------
               var(Residual) |    .134201   .0048163      .1250856    .1439806
-----------------------------+------------------------------------------------
LR test vs. linear model: chi2(2) = 203.76                Prob > chi2 = 0.0000


Note: LR test is conservative and provided only for reference.


.
. * pig data, binary outcome
. use pig_adg, clear
(Growth Performance and Abattoir Surveillance Data on Pigs, Real Data)

. generate ar_g1=ar>1 if ar~=.

. * ordinary logistic regression
. logit pn ar_g1

Iteration 0:   log likelihood = -234.95215
Iteration 1:   log likelihood = -230.59287
Iteration 2:   log likelihood = -230.59173
Iteration 3:   log likelihood = -230.59173

Logistic regression                                     Number of obs    = 341
                                                        LR chi2(1)    = 8.72
                                                        Prob > chi2    = 0.0031
Log likelihood = -230.59173                             Pseudo R2     = 0.0186

-------------+----------------------------------------------------------------
           pn | Coefficient   Std. err.       z    P>|z|      [95% con f. interval]
-------------+----------------------------------------------------------------
       ar_g1 |    .6465241    .2203379      2.93    0.003      .2146697     1.078378
       _cons |   -.1448309    .1556373     -0.93    0.352     -.4498744     .1602125
-------------+----------------------------------------------------------------

. * random effects model
. melogit pn ar_g1 || farm:

Fitting fixed-effects model:

Iteration 0:   log likelihood = -230.79889
Iteration 1:   log likelihood = -230.59177
Iteration 2:   log likelihood = -230.59173
Iteration 3:   log likelihood = -230.59173

Refining starting values:

Grid node 0:   log likelihood = -213.91338

Fitting full model:

Iteration 0:   log likelihood = -213.91338
Iteration 1:   log likelihood =    -213.52
Iteration 2:   log likelihood = -213.51177
Iteration 3:   log likelihood = -213.51176

Mixed-effects logistic regression               Number of obs     =        341
Group variable: farm                            Number of groups  =         15

                                                Obs per group:
                                                              min =         14
                                                              avg =       22.7
                                                              max =         28

Integration method: mvaghermite                 Integration pts.  =          7

                                                Wald chi2(1)      =       2.86
Log likelihood = -213.51176                     Prob > chi2       =     0.0905
-------------+----------------------------------------------------------------
           pn | Coefficient   Std. err.       z    P>|z|      [95% con f. interval]
-------------+----------------------------------------------------------------
       ar_g1 |    .4369322    .2581451      1.69    0.091     -.0690229     .9428872
       _cons |    .0196548    .3009262      0.07    0.948     -.5701497     .6094593
-------------+----------------------------------------------------------------
 farm         |
   var(_cons)|    .8774246    .4325507                     .333877     2.305861
-------------+----------------------------------------------------------------
LR test vs. logistic model: chibar2(01) = 34.16       Prob >= chibar2 =  0.0000

.
. * Reunion Island data, 3-level model for binary outcome fscr
. use reu_cfs, clear

. melogit fscr i.heifer i.ai lncfs || herd: || cow:

Fitting fixed-effects model:

Iteration 0:   log likelihood = -2011.8253
Iteration 1:   log likelihood = -2009.1421
Iteration 2:   log likelihood = -2009.1412
Iteration 3:   log likelihood = -2009.1412

Refining starting values:

Grid node 0:   log likelihood =  -2028.472

Fitting full model:

Iteration 0:   log likelihood =  -2028.472  (not concave)
Iteration 1:   log likelihood = -2016.2289  (not concave)
Iteration 2:   log likelihood = -2008.7841
Iteration 3:   log likelihood = -2000.2338
Iteration 4:   log likelihood = -1998.6784
Iteration 5:   log likelihood = -1998.6736
Iteration 6:   log likelihood = -1998.6736

Mixed-effects logistic regression               Number of obs     =      3,027

        Grouping information
        ----------------+--------------------------------------------
                        |     No. of       Observations per group
         Group variable |     groups    Minimum    Average    Maximum
        ----------------+--------------------------------------------
                   herd |         50         13       60.5        226
                    cow |      1,575          1        1.9          5
        ----------------+--------------------------------------------

Integration method: mvaghermite                 Integration pts.  =          7

                                                Wald chi2(3)      =      82.69
Log likelihood = -1998.6736                     Prob > chi2       =     0.0000
-------------+----------------------------------------------------------------
         fscr | Coefficient   Std. err.       z    P>|z|      [95% con f. interval]
-------------+----------------------------------------------------------------
      heifer |
primiparous  |   -.0620822    .0974796     -0.64    0.524     -.2531386     .1289743
             |
          ai |
         ai  |   -1.023486    .1309939     -7.81    0.000     -1.280229    -.7667427
       lncfs |    .4971351    .1019534      4.88    0.000      .2973101     .6969601
       _cons |   -1.515241    .4460314     -3.40    0.001     -2.389447    -.6410359
-------------+----------------------------------------------------------------
 herd         |
   var(_cons)|    .0732867    .0359684                    .0280068     .1917725
-------------+----------------------------------------------------------------
 herd>cow     |
   var(_cons)|    .2799605    .1219134                    .1192414     .6573041
-------------+----------------------------------------------------------------
LR test vs. logistic model: chi2(2) = 20.94               Prob > chi2 = 0.0000


Note: LR test is conservative and provided only for reference.


.
end of do-file

. exit, clear