Re: st: odd simulation results for ICC with binary outcomes

[Austin Nichols <austinnichols@gmail.com>]

Rebecca Pope <rebecca.a.pope@gmail.com>:
I don't see much difference across estimators below.
You might also want to read
 http://www.jstor.org/stable/1403259
for more options.

set seed 72114
set more off
clear all
program define betabinsim, rclass
 version 12
 syntax [, a(real 5) b(real 5) groups(int 135) n(real 100)]
 drop _all
 set obs `groups'
 gen int group = _n
 gen y1 = rbinomial(`n', rbeta(`a',`b'))
 gen y0 = `n'-y1
 reshape long y, i(group) j(event)
 xtnbreg y event, i(group) fe nolog
 scalar alph=exp(_b[_cons]+_b[event])
 scalar beta=exp(_b[_cons])
 return scalar n = `n'
 return scalar a = alph
 return scalar b = beta
 return scalar mu = alph/(alph+beta)
 return scalar icc = 1/(alph+beta+1)
 expandcl y, cluster(group event) generate(foo)
 loneway event group
 return scalar icc_anova = r(rho)
 egen double mn=mean(event), by(group)
 gen double se2=mn*(1-mn)/`n'
 su se2, mean
 scalar corrct=r(mean)
 su mn
 scalar Var=r(Var)-corrct
 ret scalar iccm=Var/(r(mean)*(1-r(mean)))
 eret clear
end
simul, r(1000): betabinsim
local truerho = 1/(1+5+5)
ttest icc==icc_anova, unpaired
ttest icc ==`truerho'
ttest icc_anova == `truerho'
ttest iccm == `truerho'
la var icc_an "ICC from ANOVA"
la var icc "ICC from xtnbreg"
la var iccm "ICC from MM"
tw kdensity icc_an||kdensity iccm||kdensity icc, xli(`=1/11')
su


On Fri, Feb 15, 2013 at 12:19 PM, Rebecca Pope <rebecca.a.pope@gmail.com> wrote:
> 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/