Notice: On March 31, it was **announced** that Statalist is moving from an email list to a **forum**. The old list will shut down on April 23, and its replacement, **statalist.org** is already up and running.

[Date Prev][Date Next][Thread Prev][Thread Next][Date Index][Thread Index]

From |
Austin Nichols <austinnichols@gmail.com> |

To |
statalist@hsphsun2.harvard.edu |

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

Date |
Fri, 15 Feb 2013 15:43:06 -0500 |

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/

**Follow-Ups**:**Re: st: odd simulation results for ICC with binary outcomes***From:*Rebecca Pope <rebecca.a.pope@gmail.com>

**References**:**st: odd simulation results for ICC with binary outcomes***From:*Rebecca Pope <rebecca.a.pope@gmail.com>

- Prev by Date:
**[no subject]** - Next by Date:
**st: RE: Graph bar shadding question.** - Previous by thread:
**st: odd simulation results for ICC with binary outcomes** - Next by thread:
**Re: st: odd simulation results for ICC with binary outcomes** - Index(es):