st: odd simulation results for ICC with binary outcomes

[Rebecca Pope <rebecca.a.pope@gmail.com>]

For background, I'm in the midst of a project to assess the
reliability of measures of physician quality. These quality measures
are scored 0 (recommended care not provided) or 1 (recommended care
provided). The powers that be suggested using -loneway- to find the
reliability of the different measures.

In the interest of full disclosure, ANOVA is not in my usual toolkit,
but my understanding is that it is a linear random effects model and
something didn't seem right about this for binary data. With
performance here being a series of binomial trials (yes/no recommended
care), the beta-binomial is a natural fit. I decided I'd run some
simulations and see how the results differ.

I focused on the intraclass correlation (ICC) for now. My thinking was
that the ICC would make the best comparison since expected values can
be calculated directly from beta shape parameters (Guimaraes, 2005)
and I would have a known true value to test against. However, while
the two estimates of the ICC are different, the ICC estimate produced
by the beta-binomial simulations is not the correct ICC while the
ANOVA analysis does converge to the correct value. The data was
generated by a beta-binomial process, so if anything should produce
that value, it should be a beta-binomial model, no?

My code is below. I ran this for different sample sizes (for which I
can also send the code if that is helpful), but except for very small
sample sizes, the results regarding the ICC are invariant to the
number of observations per group.

I would be inclined to attribute this to an as yet undiscovered error
in my code, but the a, b, mean, and variance parameters are
appropriately recovered. The mean of the beta(5,5) is 0.5. The
variance is approximately 0.0227. The ICC of a sample drawn from a
beta(5,5) should be about 0.091 (see Guimaraes, 2005). Are these
results surprising to anyone else? Does any one see any error in my
code to which this could be attributed?

Thanks so much for any help,
Rebecca

*** begin code ***
set seed 72114
set more off
capture program drop betabinsim
program define betabinsim, rclass
 version 12
 syntax [, a(real 5) b(real 5) groups(integer 135) n(real 100)]
 drop _all
 tempvar group foo event
 set obs `groups'
 gen `group' = _n // pseudo-identifiers for physician groups
 gen y1 = rbinomial(`n', rbeta(`a',`b'))  // draw successes for each
group around beta(a,b)
 gen y0 = `n'-y1  // failues for each PFSA
  // Set-up for beta-binomial
 reshape long y, i(`group') j(`event')
  // Guimaraes (2005) procedure for estimating beta-binomial parameters
 xtnbreg y `event', i(`group') fe nolog
 local alpha exp(_b[_cons]+_b[`event'])
 local beta exp(_b[_cons])

 return scalar n = `n'
 return scalar a = `alpha'
 return scalar b = `beta'
 return scalar mu = `alpha'/(`alpha'+`beta')
 return scalar icc = 1/(`alpha'+`beta'+1)

  // Set-up for ANOVA
 expandcl y, cluster(`group' `event') generate(`foo')
  // Run ANOVA & save ICC
 loneway `event' `group'
 return scalar icc_anova = `r(rho)'
end
simulate sampsi=r(n) alpha=r(a) beta=r(b) mu=r(mu) icc=r(icc)
icc_anova=r(icc_anova), reps(10000): betabinsim

local truerho = 1/(1+5+5)
ttest icc==icc_anova, unpaired
ttest icc ==`truerho'
ttest icc_anova == `truerho'
*** end ***

References:
Guimaraes, P. (2005) A simple approach to fit the beta-binomial model.
The Stata Journal, 5(3), 385-394. Available at:
http://www.stata-journal.com/article.html?article=st0089
*
*   For searches and help try:
*   http://www.stata.com/help.cgi?search
*   http://www.stata.com/support/faqs/resources/statalist-faq/
*   http://www.ats.ucla.edu/stat/stata/