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

. * do-file for part of lecture 11a of VHM 802/812, Winter 2015
. version 13

. set more off

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

.
. * 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      =      2178
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 |      Coef.    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 regression: 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      =      2178
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 |      Coef.    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 regression: 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

. 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      =      2178
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 |      Coef.    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 regression: 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.8919
Iteration 2:   log restricted-likelihood = -1495.8804
Iteration 3:   log restricted-likelihood = -1495.8804

Computing standard errors:

Mixed-effects REML regression                   Number of obs      =      3027

----------------+------------------------------------------
                |   No. of       Observations per Group
 Group Variable |   Groups    Minimum    Average    Maximum
----------------+------------------------------------------
           herd |       50         13       60.5        226
            cow |     1575          1        1.9          5
----------------+------------------------------------------

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

-------------+----------------------------------------------------------------
        lncfs |      Coef.    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      .0134606    .0298868
-----------------------------+------------------------------------------------
               var(Residual) |   .1341146   .0048118      .1250077     .143885
-----------------------------+------------------------------------------------
LR test vs. linear regression:       chi2(2) =   221.39   Prob > chi2 = 0.0000


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


. * 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     .1546348     .2665226
 -----------------------------+------------------------------------------------

. * model checking
. predict res_lact, rstandard

. summarize res_lact, d

                   Standardized residuals
-------------------------------------------------------------
      Percentiles      Smallest
 1%    -1.902202      -2.385942
 5%    -1.376046      -2.376865
10%    -1.118899       -2.24191       Obs                3027
25%    -.6474585      -2.235393       Sum of Wgt.        3027

50%    -.0968802                      Mean          -1.91e-10
                        Largest       Std. Dev.      .9372907
75%     .5640011       3.537244
90%     1.182198       3.713121       Variance       .8785139
95%     1.670863       4.115439       Skewness       .6035033
99%     2.719746       4.525204       Kurtosis       3.792215

. qnorm res_lact

. histogram res_lact
(bin=34, start=-2.3859422, width=.20326901)

. swilk res_lact

                   Shapiro-Wilk W test for normal data

    Variable |    Obs       W           V         z       Prob>z
-------------+--------------------------------------------------
    res_lact |   3027    0.98077     33.204     9.042    0.00000

. predict fitted, fit /* includes all random effects */

. scatter res_lact fitted

. * predicted random effects at higher levels
. predict ref*, reffects

. bysort cow: generate within_c=_n

. bysort herd: generate within_h=_n

. * herd-level random effects
. summarize ref1 if within_h==1, d

                  BLUP r.e. for herd: _cons
-------------------------------------------------------------
      Percentiles      Smallest
 1%    -.2785034      -.2785034
 5%    -.1536851       -.183448
10%    -.1307647      -.1536851       Obs                  50
25%     -.079566      -.1500316       Sum of Wgt.          50

50%    -.0008414                      Mean           1.12e-10
                        Largest       Std. Dev.      .1074885
75%     .0726699       .1262656
90%     .1181743       .1290758       Variance       .0115538
95%     .1290758       .2022623       Skewness       .2345923
99%     .3407647       .3407647       Kurtosis       4.076068

. qnorm ref1 if within_h==1

. swilk ref1 if within_h==1

                   Shapiro-Wilk W test for normal data

    Variable |    Obs       W           V         z       Prob>z
-------------+--------------------------------------------------
        ref1 |     50    0.97805      1.032     0.068    0.47286

. * cow-level random effects
. summarize ref2 if within_c==1, d

                  BLUP r.e. for cow: _cons
-------------------------------------------------------------
      Percentiles      Smallest
 1%    -.1376604      -.1763716
 5%    -.1015879      -.1629009
10%     -.078711       -.160541       Obs                1575
25%    -.0428026      -.1563797       Sum of Wgt.        1575

50%    -.0046796                      Mean           4.10e-11
                        Largest       Std. Dev.      .0646133
75%     .0381531       .2430278
90%     .0797535       .2502964       Variance       .0041749
95%     .1140806       .2507932       Skewness       .4459166
99%     .1757991       .2671055       Kurtosis       3.654445

. qnorm ref2 if within_c==1

. swilk ref2 if within_c==1

                   Shapiro-Wilk W test for normal data

    Variable |    Obs       W           V         z       Prob>z
-------------+--------------------------------------------------
        ref2 |   1575    0.98805     11.402     6.136    0.00000

.
end of do-file

. exit, clear