-gllamm- and ICC

[Jen Ward]

Hi,

I am analysing data from a cluster randomised trial using -gllamm- and to estimate relative risks I am specifying binomial family with log link; I would like to estimate the ICC but I am not sure if/how this can be done - can anyone advise if there is a formula?

Thanks,
Jen

[Andrew Musau]

gllamm is from SSC, as you are asked to explain in FAQ Advice #12. For a 3-level mixed-effects logit model, you can calculate ICC as follows:

$$\text{ICC}_{1}= \frac{Var(\_cons_{1})}{Var(\_cons_{1})+ Var(\_cons_{2}) + \frac{\text{pi}^{2}}{3}}$$

and

$$\text{ICC}_{2}=\text{ICC}_{1}+ \frac{Var(\_cons_{2})}{Var(\_cons_{1})+ Var(\_cons_{2}) + \frac{\text{pi}^{2}}{3}}$$

The gllamm coefficients corresponding to the variance components stored in r(table) need to be squared to obtain the variances of the random effects. Here is a model estimated using melogit where ICCs are calculated using estat icc and the same model estimated using gllamm and the ICCs computed by hand.

webuse tvsfpors, clear
gen thkbin= thk>2
melogit thkbin prethk cc##tv || school: || class:, nolog
estat icc

gen cctv= c.cc#c.tv
gllamm thkbin prethk cc tv cctv, i(class school) fam(binomial) link(logit) nolog
scalar icc1= r(table)["b", "scho2:_cons"]^2/(r(table)["b", "clas1:_cons"]^2+  r(table)["b", "scho2:_cons"]^2 + (c(pi)^2)/3)
di icc1

scalar icc2= icc1+(r(table)["b", "clas1:_cons"]^2/(r(table)["b", "clas1:_cons"]^2+  r(table)["b", "scho2:_cons"]^2 + (c(pi)^2)/3))
di icc2

Res.:

. melogit thkbin prethk cc##tv || school: || class:, nolog

Mixed-effects logistic regression               Number of obs     =      1,600

        Grouping information
        -------------------------------------------------------------
                        |     No. of       Observations per group
         Group variable |     groups    Minimum    Average    Maximum
        ----------------+--------------------------------------------
                 school |         28         18       57.1        137
                  class |        135          1       11.9         28
        -------------------------------------------------------------

Integration method: mvaghermite                 Integration pts.  =          7

                                                Wald chi2(4)      =      91.42
Log likelihood = -1027.851                      Prob > chi2       =     0.0000
------------------------------------------------------------------------------
      thkbin | Coefficient  Std. err.      z    P>|z|     [95% conf. interval]
-------------+----------------------------------------------------------------
      prethk |    .395364   .0463324     8.53   0.000     .3045541    .4861738
        1.cc |   1.038282   .2447573     4.24   0.000     .5585667    1.517998
        1.tv |   .3325112    .235739     1.41   0.158    -.1295289    .7945512
             |
       cc#tv |
        1 1  |  -.4643693   .3426676    -1.36   0.175    -1.135986    .2072469
             |
       _cons |  -1.246478   .1956927    -6.37   0.000    -1.630029   -.8629273
-------------+----------------------------------------------------------------
school       |
   var(_cons)|   .0628837   .0616755                       .009198    .4299174
-------------+----------------------------------------------------------------
school>class |
   var(_cons)|   .1649057   .0813311                      .0627227    .4335572
------------------------------------------------------------------------------
LR test vs. logistic model: chi2(2) = 17.61               Prob > chi2 = 0.0001

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

. 
. estat icc

Residual intraclass correlation

------------------------------------------------------------------------------
                       Level |        ICC   Std. err.     [95% conf. interval]
-----------------------------+------------------------------------------------
                      school |   .0178766   .0173592      .0026144    .1122118
                class|school |    .064756   .0224644      .0323835    .1252994
------------------------------------------------------------------------------

. 
. 
. 
. gen cctv= c.cc#c.tv

. 
. gllamm thkbin prethk cc tv cctv, i(class school) fam(binomial) link(logit) nolog
 
number of level 1 units = 1600
number of level 2 units = 135
number of level 3 units = 28
 
Condition Number = 15.417248
 
gllamm model 
 
log likelihood = -1027.8514
 
------------------------------------------------------------------------------
      thkbin | Coefficient  Std. err.      z    P>|z|     [95% conf. interval]
-------------+----------------------------------------------------------------
      prethk |   .3953702   .0463311     8.53   0.000     .3045629    .4861775
          cc |   1.038113   .2447844     4.24   0.000     .5583449    1.517882
          tv |   .3322967   .2358046     1.41   0.159    -.1298718    .7944653
        cctv |  -.4642312   .3426383    -1.35   0.175     -1.13579    .2073275
       _cons |  -1.246279   .1959921    -6.36   0.000    -1.630417   -.8621417
------------------------------------------------------------------------------
 
 
Variances and covariances of random effects
------------------------------------------------------------------------------

 
***level 2 (class)
 
    var(1): .16496417 (.0813068)
 
***level 3 (school)
 
    var(1): .06279013 (.06147452)
------------------------------------------------------------------------------

 

. 
. scalar icc1= r(table)["b", "scho2:_cons"]^2/(r(table)["b", "clas1:_cons"]^2+  r(table)["b", "scho2:_cons"]^2 + (c(pi)^2)
> /3)

. 
. di icc1
.01785016

. 
. 
. 
. scalar icc2= icc1+(r(table)["b", "clas1:_cons"]^2/(r(table)["b", "clas1:_cons"]^2+  r(table)["b", "scho2:_cons"]^2 + (c(
> pi)^2)/3))

. 
. di icc2
.06474666

. 

[Leonardo Guizzetti]

Andrew can this be done using the binomial family with log link, as the OP requested? I'm not certain that it would be the same, as the (pi^2)/3 constant is the variance of he logistic distribution.

[Andrew Musau]

Correct, #2 only applies to -link(logit)-. I somehow missed the -link(log)- part. In the linear case where ICC= var(_cons_j)/total variance and with a logit link, the computation is straight-forward. I have no idea of how one proceeds with a log link.