Estimating risk difference for binary outcome in cluster RCT

[Sonia Lewis]

I am analysing data from a cluster RCT looking at the effect of an intervention (int) in health centres (centre) to reduce inappropriate antibiotic prescribing (abu). We initially planned to present results as ORs using mixed models with random effects for health centre. This is the code:

meglm abu int || centre: , family(binomial) link(logit) eform

However, antibiotic prescribing is over 90% and the OR showed a very strong effect (OR=0.25) while the absolute difference was small (5%). This seems misleading, and we thought maybe this was because the outcome is not rare, so we decided to present the results as RRs rather than ORs using this code:

meglm abu int || centre: , family(poisson) link(log) eform

This gives an RR of 0.93, which seems more reasonable.

We also wanted to present Risk Difference (RD) and have tried a few different ways to get this. Binomial with identity or log link as suggested in this post

https://www.statalist.org/forums/forum/general-stata-discussion/general/1693010-adjusted-risk-ratios-and-risk-differences-after-logistic-regression

gives error messages saying that identity and log links are not allowed with binomial distribution. We also tried:

meglm abu int || centre: , family(gaussian) link(identity)

But this gives an error saying that initial values are not feasible.

Can anyone suggest another way of doing this? Or is the best way to use margins?

[Joseph Coveney]

We also wanted to present Risk Difference (RD) . . . But this gives an error saying that initial values are not feasible.

Can anyone suggest another way of doing this? Or is the best way to use margins?

-margins- is probably your best bet, but you can fit a linear model to your data using a noniterative method, such as -regress- or -xtreg , re-, and then use -nlcom- to obtain all three treatment-effect measures, odds ratio (OR), risk ratio (RR) as well as RD. I show examples below of both (begin at the "Begin here" comment; the top part is just to create a toy dataset for use in illustration), along with an approach to obtain all three with a random-effects logistic regression model via -margins- (and -nlcom- for RR). They all give about the same respective values for each of the three treatment-effect measures.

The fact that -meglm- cannot find feasible initial values is a bit troubling in that it suggests that you have relatively few clusters, or that they're relatively small, and at least a significant fraction of them have all the same outcome: all members of the cluster—presumably assigned to the control-treatment condition—are positive for inappropriate prescribing practices.

.ÿ
.ÿversionÿ17.0

.ÿ
.ÿclearÿ*

.ÿ
.ÿ//ÿseedem
.ÿsetÿseedÿ976106965

.ÿ
.ÿquietlyÿsetÿobsÿ100

.ÿgenerateÿbyteÿsidÿ=ÿ_n

.ÿgenerateÿdoubleÿsid_uÿ=ÿrnormal()

.ÿ
.ÿgenerateÿbyteÿtrtÿ=ÿmod(_n,ÿ2)

.ÿgenerateÿdoubleÿxbuÿ=ÿlogit(0.982ÿ-ÿtrtÿ*ÿ0.05)ÿ+ÿsid_u

.ÿ
.ÿquietlyÿexpandÿ100

.ÿgenerateÿbyteÿabuÿ=ÿrbinomial(1,ÿinvlogit(xbu))

.ÿ
.ÿ*
.ÿ*ÿBeginÿhere
.ÿ*
.ÿxtlogitÿabuÿi.trt,ÿi(sid)ÿreÿorÿnolog

Random-effectsÿlogisticÿregressionÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿNumberÿofÿobsÿÿÿÿ=ÿ10,000
Groupÿvariable:ÿsidÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿNumberÿofÿgroupsÿ=ÿÿÿÿ100

Randomÿeffectsÿu_iÿ~ÿGaussianÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿObsÿperÿgroup:
ÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿminÿ=ÿÿÿÿ100
ÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿavgÿ=ÿÿ100.0
ÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿmaxÿ=ÿÿÿÿ100

Integrationÿmethod:ÿmvaghermiteÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿIntegrationÿpts.ÿ=ÿÿÿÿÿ12

ÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿWaldÿchi2(1)ÿÿÿÿÿ=ÿÿ36.82
Logÿlikelihoodÿ=ÿ-2284.7862ÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿProbÿ>ÿchi2ÿÿÿÿÿÿ=ÿ0.0000

------------------------------------------------------------------------------
ÿÿÿÿÿÿÿÿÿabuÿ|ÿOddsÿratioÿÿÿStd.ÿerr.ÿÿÿÿÿÿzÿÿÿÿP>|z|ÿÿÿÿÿ[95%ÿconf.ÿinterval]
-------------+----------------------------------------------------------------
ÿÿÿÿÿÿÿ1.trtÿ|ÿÿÿ.2507516ÿÿÿ.0571619ÿÿÿÿ-6.07ÿÿÿ0.000ÿÿÿÿÿ.1603991ÿÿÿÿ.3919997
ÿÿÿÿÿÿÿ_consÿ|ÿÿÿ45.13468ÿÿÿ7.943071ÿÿÿÿ21.65ÿÿÿ0.000ÿÿÿÿÿ31.96766ÿÿÿÿ63.72501
-------------+----------------------------------------------------------------
ÿÿÿÿ/lnsig2uÿ|ÿÿ-.0507031ÿÿÿ.2022377ÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿ-.4470816ÿÿÿÿ.3456754
-------------+----------------------------------------------------------------
ÿÿÿÿÿsigma_uÿ|ÿÿÿ.9749671ÿÿÿ.0985875ÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿ.7996823ÿÿÿÿ1.188673
ÿÿÿÿÿÿÿÿÿrhoÿ|ÿÿÿ.2241662ÿÿÿ.0351723ÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿ.162747ÿÿÿÿ.3004466
------------------------------------------------------------------------------
Note:ÿEstimatesÿareÿtransformedÿonlyÿinÿtheÿfirstÿequationÿtoÿoddsÿratios.
Note:ÿ_consÿestimatesÿbaselineÿoddsÿ(conditionalÿonÿzeroÿrandomÿeffects).
LRÿtestÿofÿrho=0:ÿchibar2(01)ÿ=ÿ326.46ÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿProbÿ>=ÿchibar2ÿ=ÿ0.000

.ÿ//ÿHere:
.ÿmarginsÿ,ÿdydx(trt)

ConditionalÿmarginalÿeffectsÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿNumberÿofÿobsÿ=ÿ10,000
ModelÿVCE:ÿOIM

Expression:ÿPr(abu=1),ÿpredict(pr)
dy/dxÿwrt:ÿÿ1.trt

------------------------------------------------------------------------------
ÿÿÿÿÿÿÿÿÿÿÿÿÿ|ÿÿÿÿÿÿÿÿÿÿÿÿDelta-method
ÿÿÿÿÿÿÿÿÿÿÿÿÿ|ÿÿÿÿÿÿdy/dxÿÿÿstd.ÿerr.ÿÿÿÿÿÿzÿÿÿÿP>|z|ÿÿÿÿÿ[95%ÿconf.ÿinterval]
-------------+----------------------------------------------------------------
ÿÿÿÿÿÿÿ1.trtÿ|ÿÿ-.0774658ÿÿÿ.0145561ÿÿÿÿ-5.32ÿÿÿ0.000ÿÿÿÿ-.1059951ÿÿÿ-.0489364
------------------------------------------------------------------------------
Note:ÿdy/dxÿforÿfactorÿlevelsÿisÿtheÿdiscreteÿchangeÿfromÿtheÿbaseÿlevel.

.ÿ
.ÿ//ÿAlso,ÿthere'sÿthis:
.ÿquietlyÿmarginsÿtrt,ÿpost

.ÿnlcomÿ(RR:ÿ_b[1.trt]ÿ/ÿ_b[0.trt])ÿ(RD:ÿ_b[1.trt]ÿ-ÿ_b[0.trt]),ÿnoheader

------------------------------------------------------------------------------
ÿÿÿÿÿÿÿÿÿÿÿÿÿ|ÿCoefficientÿÿStd.ÿerr.ÿÿÿÿÿÿzÿÿÿÿP>|z|ÿÿÿÿÿ[95%ÿconf.ÿinterval]
-------------+----------------------------------------------------------------
ÿÿÿÿÿÿÿÿÿÿRRÿ|ÿÿÿ.9198995ÿÿÿ.0149219ÿÿÿÿ61.65ÿÿÿ0.000ÿÿÿÿÿ.8906532ÿÿÿÿ.9491459
ÿÿÿÿÿÿÿÿÿÿRDÿ|ÿÿ-.0774658ÿÿÿ.0145561ÿÿÿÿ-5.32ÿÿÿ0.000ÿÿÿÿ-.1059951ÿÿÿ-.0489364
------------------------------------------------------------------------------

.ÿ
.ÿ//ÿAndÿthis:
.ÿregressÿabuÿibn.trt,ÿnoconstantÿvce(clusterÿsid)

LinearÿregressionÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿNumberÿofÿobsÿÿÿÿÿ=ÿÿÿÿÿ10,000
ÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿF(2,ÿ99)ÿÿÿÿÿÿÿÿÿÿ=ÿÿÿ20274.30
ÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿProbÿ>ÿFÿÿÿÿÿÿÿÿÿÿ=ÿÿÿÿÿ0.0000
ÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿR-squaredÿÿÿÿÿÿÿÿÿ=ÿÿÿÿÿ0.9305
ÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿRootÿMSEÿÿÿÿÿÿÿÿÿÿ=ÿÿÿÿÿ.25415

ÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿ(Std.ÿerr.ÿadjustedÿforÿ100ÿclustersÿinÿsid)
------------------------------------------------------------------------------
ÿÿÿÿÿÿÿÿÿÿÿÿÿ|ÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿRobust
ÿÿÿÿÿÿÿÿÿabuÿ|ÿCoefficientÿÿstd.ÿerr.ÿÿÿÿÿÿtÿÿÿÿP>|t|ÿÿÿÿÿ[95%ÿconf.ÿinterval]
-------------+----------------------------------------------------------------
ÿÿÿÿÿÿÿÿÿtrtÿ|
ÿÿÿÿÿÿÿÿÿÿ0ÿÿ|ÿÿÿÿÿÿ.9672ÿÿÿ.0051174ÿÿÿ189.00ÿÿÿ0.000ÿÿÿÿÿÿ.957046ÿÿÿÿÿ.977354
ÿÿÿÿÿÿÿÿÿÿ1ÿÿ|ÿÿÿÿÿÿ.8906ÿÿÿ.0128192ÿÿÿÿ69.47ÿÿÿ0.000ÿÿÿÿÿ.8651639ÿÿÿÿ.9160361
------------------------------------------------------------------------------

.ÿnlcomÿ///
>ÿÿÿÿÿÿÿÿÿ(OR:ÿexp(logit(_b[1.trt])ÿ-ÿlogit(_b[0.trt])))ÿ///
>ÿÿÿÿÿÿÿÿÿ(RR:ÿ_b[1.trt]ÿ/ÿ_b[0.trt])ÿ(RD:ÿ_b[1.trt]ÿ-ÿ_b[0.trt]),ÿnoheader

------------------------------------------------------------------------------
ÿÿÿÿÿÿÿÿÿabuÿ|ÿCoefficientÿÿStd.ÿerr.ÿÿÿÿÿÿzÿÿÿÿP>|z|ÿÿÿÿÿ[95%ÿconf.ÿinterval]
-------------+----------------------------------------------------------------
ÿÿÿÿÿÿÿÿÿÿORÿ|ÿÿÿ.2760724ÿÿÿ.0574679ÿÿÿÿÿ4.80ÿÿÿ0.000ÿÿÿÿÿ.1634374ÿÿÿÿ.3887074
ÿÿÿÿÿÿÿÿÿÿRRÿ|ÿÿÿ.9208023ÿÿÿÿ.014121ÿÿÿÿ65.21ÿÿÿ0.000ÿÿÿÿÿ.8931257ÿÿÿÿÿ.948479
ÿÿÿÿÿÿÿÿÿÿRDÿ|ÿÿÿÿÿ-.0766ÿÿÿ.0138029ÿÿÿÿ-5.55ÿÿÿ0.000ÿÿÿÿ-.1036532ÿÿÿ-.0495468
------------------------------------------------------------------------------

.ÿ
.ÿ//ÿOrÿthis:
.ÿxtregÿabuÿi.trt,ÿi(sid)ÿreÿvce(robust)

Random-effectsÿGLSÿregressionÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿNumberÿofÿobsÿÿÿÿÿ=ÿÿÿÿÿ10,000
Groupÿvariable:ÿsidÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿNumberÿofÿgroupsÿÿ=ÿÿÿÿÿÿÿÿ100

R-squared:ÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿObsÿperÿgroup:
ÿÿÿÿÿWithinÿÿ=ÿ0.0000ÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿminÿ=ÿÿÿÿÿÿÿÿ100
ÿÿÿÿÿBetweenÿ=ÿ0.2373ÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿavgÿ=ÿÿÿÿÿÿ100.0
ÿÿÿÿÿOverallÿ=ÿ0.0222ÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿmaxÿ=ÿÿÿÿÿÿÿÿ100

ÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿWaldÿchi2(1)ÿÿÿÿÿÿ=ÿÿÿÿÿÿ30.80
corr(u_i,ÿX)ÿ=ÿ0ÿ(assumed)ÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿProbÿ>ÿchi2ÿÿÿÿÿÿÿ=ÿÿÿÿÿ0.0000

ÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿ(Std.ÿerr.ÿadjustedÿforÿ100ÿclustersÿinÿsid)
------------------------------------------------------------------------------
ÿÿÿÿÿÿÿÿÿÿÿÿÿ|ÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿRobust
ÿÿÿÿÿÿÿÿÿabuÿ|ÿCoefficientÿÿstd.ÿerr.ÿÿÿÿÿÿzÿÿÿÿP>|z|ÿÿÿÿÿ[95%ÿconf.ÿinterval]
-------------+----------------------------------------------------------------
ÿÿÿÿÿÿÿ1.trtÿ|ÿÿÿÿÿ-.0766ÿÿÿ.0138029ÿÿÿÿ-5.55ÿÿÿ0.000ÿÿÿÿ-.1036532ÿÿÿ-.0495468
ÿÿÿÿÿÿÿ_consÿ|ÿÿÿÿÿÿ.9672ÿÿÿ.0051174ÿÿÿ189.00ÿÿÿ0.000ÿÿÿÿÿ.9571701ÿÿÿÿ.9772299
-------------+----------------------------------------------------------------
ÿÿÿÿÿsigma_uÿ|ÿÿ.06485712
ÿÿÿÿÿsigma_eÿ|ÿÿ.24590176
ÿÿÿÿÿÿÿÿÿrhoÿ|ÿÿ.06504063ÿÿÿ(fractionÿofÿvarianceÿdueÿtoÿu_i)
------------------------------------------------------------------------------

.ÿnlcomÿ///
>ÿÿÿÿÿÿÿÿÿ(OR:ÿexp(logit(_b[1.trt]ÿ+ÿ_b[_cons])ÿ-ÿlogit(_b[_cons])))ÿ///
>ÿÿÿÿÿÿÿÿÿ(RR:ÿ(_b[1.trt]ÿ+ÿ_b[_cons])ÿ/ÿ_b[_cons]),ÿnoheader

------------------------------------------------------------------------------
ÿÿÿÿÿÿÿÿÿabuÿ|ÿCoefficientÿÿStd.ÿerr.ÿÿÿÿÿÿzÿÿÿÿP>|z|ÿÿÿÿÿ[95%ÿconf.ÿinterval]
-------------+----------------------------------------------------------------
ÿÿÿÿÿÿÿÿÿÿORÿ|ÿÿÿ.2760724ÿÿÿ.0574679ÿÿÿÿÿ4.80ÿÿÿ0.000ÿÿÿÿÿ.1634374ÿÿÿÿ.3887074
ÿÿÿÿÿÿÿÿÿÿRRÿ|ÿÿÿ.9208023ÿÿÿÿ.014121ÿÿÿÿ65.21ÿÿÿ0.000ÿÿÿÿÿ.8931257ÿÿÿÿÿ.948479
------------------------------------------------------------------------------

.ÿ
.ÿexit

endÿofÿdo-file


.

Summary of fitted coefficients of treatment effect:

ÿÿÿMixed-effectsÿLogisticÿÿLinearÿRegressionÿÿRandom-effectsÿLinearÿÿÿÿÿÿÿÿÿÿÿGEE*
ORÿÿÿÿÿÿ0.25ÿ±ÿ0.06ÿÿÿÿÿÿÿÿÿÿÿ0.27ÿ±ÿ0.06ÿÿÿÿÿÿÿÿÿÿÿ0.28ÿ±ÿ0.06ÿÿÿÿÿÿÿÿÿÿÿ0.28ÿ±ÿ0.06
RRÿÿÿÿÿÿ0.92ÿ±ÿ0.01ÿÿÿÿÿÿÿÿÿÿÿ0.92ÿ±ÿ0.01ÿÿÿÿÿÿÿÿÿÿÿ0.92ÿ±ÿ0.01ÿÿÿÿÿÿÿÿÿÿÿ0.92ÿ±ÿ0.01
RDÿÿÿÿ-0.077ÿ±ÿ0.015ÿÿÿÿÿÿÿÿ-0.077ÿ±ÿ0.014ÿÿÿÿÿÿÿÿ-0.077ÿ±ÿ0.014ÿÿÿÿÿÿÿÿ-0.077ÿ±ÿ0.014

*ÿOutputÿnotÿshown,ÿbutÿobtainedÿfrom:
xtgeeÿabuÿi.trt,ÿi(sid)ÿfamily(binomial)ÿlink(logit)ÿcorr(exchangeable)ÿeformÿnolog
xtgeeÿabuÿi.trt,ÿi(sid)ÿfamily(binomial)ÿlink(log)ÿcorr(exchangeable)ÿeformÿnolog
xtgeeÿabuÿi.trt,ÿi(sid)ÿfamily(binomial)ÿlink(identity)ÿcorr(exchangeable)ÿvce(robust)ÿnolog

[Sonia Lewis]

Thanks so much Joseph, that's really helpful. I tried those commands with a preliminary version of the real dataset and also added mixed-effects Poisson model using:

xtpoisson abu i.armx, i(locationx) re irr nolog
quietly margins armx, post
nlcom (RD: _b[1.armx] - _b[0.armx]), noheader

In your example, the estimates from the different models are all similar, so how would you decide which model to use? With the real data, there are differences between the estimates depending on which model is used, and I'm not sure why, or which one is the most appropriate. To give a little more information about the dataset, there are 48 clusters with between ~100 and ~4000 patients. No clusters have 100% of patients with inappropriate antibiotic use, but 7 clusters have more than 99% (i.e. <10 patients not receiving antibiotics), 6 of which are in the control arm. Could this be why the different models show different effects? If so, which is the most appropriate one to use with this data? And should RR be reported, or is it ok to use OR?

[Joseph Coveney]

With the real data, there are differences between the estimates depending on which model is used, and I'm not sure why, or which one is the most appropriate.

It is a curious observation. No guarantee, but there is something that you could do that might help to clear it up: with your dataset, execute the following

contract int centre abu, freq(count)

and then use -dataex- to insert the 96 rows that result into a new post. (Omit any generated label define commands from the -dataex- output that might contain centre identifiers.)

Those are the variable names that are used in your first post. If the variable names are now as shown in your second, then

contract armx locationx abu, freq(count)

followed by -dataex-, which will list up to 100 rows of data without requiring any tweaking of its options.

[Sonia Lewis]

Happy new year!

Here is the data extracted using your code:

. contract armx locationx abu, freq(count)

. dataex

----------------------- copy starting from the next line -----------------------

* Example generated by -dataex-. For more info, type help dataex
clear
input byte abu long(armx locationx) int count
0 0  1    2
1 0  1  712
0 0  4   39
1 0  4  711
0 0  5   72
1 0  5  753
0 0  8   90
1 0  8 1325
0 0  9   38
1 0  9  975
0 0 10    2
1 0 10 1074
0 0 11   27
1 0 11 1263
0 0 15   13
1 0 15  420
0 0 16   29
1 0 16  609
0 0 20   11
1 0 20  249
0 0 23   62
1 0 23  529
0 0 26  506
1 0 26  367
0 0 27   10
1 0 27  475
0 0 29    1
1 0 29  861
0 0 31    2
1 0 31 1213
0 0 33  250
1 0 33  735
0 0 34   61
1 0 34  438
0 0 35   17
1 0 35  764
0 0 36   10
1 0 36  776
0 0 39  209
1 0 39 1349
0 0 40   39
1 0 40 1874
0 0 42    4
1 0 42 1037
0 0 44   27
1 0 44 1357
0 0 48    2
1 0 48  664
0 1  2   48
1 1  2   70
0 1  3   84
1 1  3  933
0 1  6    7
1 1  6  301
0 1  7   16
1 1  7  150
0 1 12  109
1 1 12  279
0 1 13    7
1 1 13  540
0 1 14   22
1 1 14 4411
0 1 17   14
1 1 17  684
0 1 18   66
1 1 18  227
0 1 19  107
1 1 19  899
0 1 21  116
1 1 21  375
0 1 22   24
1 1 22  512
0 1 24   67
1 1 24  441
0 1 25  175
1 1 25  678
0 1 28   97
1 1 28   62
0 1 30  147
1 1 30 1155
0 1 32   78
1 1 32  262
0 1 37  121
1 1 37 1184
0 1 38   91
1 1 38  686
0 1 41   94
1 1 41  631
0 1 43   71
1 1 43 1184
0 1 45   60
1 1 45  245
0 1 46   28
1 1 46  987
0 1 47  142
1 1 47 1702
end
label values armx armx
label def armx 0 "Control", modify
label def armx 1 "Intervention", modify

I have done some further exploration using the different modelling approaches with margins and nlcom you suggested, and have some more questions. The preferred option from your suggestions was mixed-effects logistic with xtlogit, as this fits best with our statistical analysis plan. However there are some problems with this:

1. This model has problems when covariates are added and reports the error "Initial values not feasible". For example when a binary variable for whether patients are under or over 16 years.

xtlogit abu i.armx i.age_cat16x districtx, i(locationx) re or nolog

I assume this is because some clusters don't have patients in all of the arm, abu, and age group combinations (I checked and this is the case - 187 combinations instead of 192). I tried running a similar model with meglm using this code:

meglm abu i.armx i.age_cat16x districtx || locationx: , family(binomial) link(logit) eform

and this model does not have the same problems when covariates. ORs are the same for xtlogit and meglm, but confidence intervals and p-values differ slightly. Is it ok to use meglm instead of xtlogit?


2. When using the RR after nlcom with both xtlogit and meglm, the p-values seem to be incorrect. For meglm, overall, the RR and RD are like this:

------------------------------------------------------------------------------
| Coefficient Std. err. z P&gt;|z| [95% conf. interval]
-------------+----------------------------------------------------------------
RR | .8889392 .0375048 23.70 0.000 .8154312 .9624472
RD | -.1044098 .0358089 -2.92 0.004 -.1745938 -.0342257
------------------------------------------------------------------------------

This p-value seems too small compared to the confidence interval. For one of the sub-group models the 95% CI overlaps 1, while the p-value is still very small:

------------------------------------------------------------------------------
| Coefficient Std. err. z P&gt;|z| [95% conf. interval]
-------------+----------------------------------------------------------------
RR | .9538965 .0371261 25.69 0.000 .8811306 1.026662
RD | -.0432947 .0352986 -1.23 0.220 -.1124787 .0258892
------------------------------------------------------------------------------

Any ideas why this might be and how to get the correct p-values for RRs? The p-values for the RDs seem to be ok.


3. When testing for interaction between arm and other variables using LRT, these are calculated for the underlying model, not the RR derived from nlcom. For mixed-effects logistic models most interactions appear to be highly significant because of the strong effects sizes with ORs, whereas the effect sizes and stratum-specific effects are much smaller when presenting RRs. So to get the RRs directly from the same models that are used to test for interactions, I tried xtpoisson

xtpoisson abu i.armx i.age_cat16x districtx, i(locationx) re irr

and also meglm with poisson and log link

meglm abu i.armx i.age_cat16x districtx || locationx: , family(poisson) link(log) eform

These both work ok and the results are similar to each other. Is it appropriate to use poisson models with these data, and if so, which is better?


4. I also tried the linear regression (regress) and random effects linear regression (xtreg) models with nlcom, but the p-values for the RRs derived through nlcom have the same problem as the RRs from logistic models in point 2. and don't run for some sub-group analyses. Also, LRT does not work with vce(robust).



Any advice on which model is most appropriate and which effect measure is most appropriate with these data would be much appreciated!

[Sonia Lewis]

For point 2 in my previous post, I calculated the p-values for the RRs using the 95% Confidence Intervals and the procedure in this paper: https://www.bmj.com/content/343/bmj.d2304 . The overall p-value in point 2 is RR 0.89 (0.82-96), p=0.0053, and the one for the sub-group is RR 0.95 (0.88-1.03), p=0.22. This looks ok. But I'm still not sure why the nlcom output gives appropriate 95% CI but wrong p-value. Any ideas?

[Joseph Coveney]

For point 2 . . . I'm still not sure why the nlcom output gives appropriate 95% CI but wrong p-value. Any ideas?

Yeah, the p-value from -nlcom- is for the default test of the null hypothesis that the computed value (RR in this case) equals zero. -nlcom- doesn't know that you want the test that RR equals one.

You can do the following.

meglm . . .
margins trt, post
nlcom (RR: _b[1.trt] / _b[0.trt]), post
test RR = 1

Thanks for posting the dataset. I do have suggestions and comments on your other questions, but it will be probably a day before I'll have time to get back to you with them.

[Joseph Coveney]

Sorry, I accidentally stuck with my earlier variable name for the intervention. For the code just above, substitute yours (armx) for mine (trt).

[Joseph Coveney]

1. This model has problems when covariates are added and reports the error "Initial values not feasible". . . . I assume this is because some clusters don't have patients in all of the arm, abu, and age group combinations

That's possible, but what is districtx? Without a factor-variable prefix, it's handled by default as a continuous predictor, but from its name it sounds as if it is categorical. If it has more than two categories coded as zero and one, then this might be misspecified.

I tried running a similar model with meglm . . . and this model does not have the same problems when covariates. ORs are the same for xtlogit and meglm, but confidence intervals and p-values differ slightly. Is it ok to use meglm instead of xtlogit?

The difference in confidence bounds and p-values that you see between xtlogit , re and meglm , family(binomial) link(logit) are likely because by their respective defaults the former uses 12 integration points and the latter only seven. (This might also be the basis for the difficulty that you experience with the former.) They should behave identically if you set the number of integration points the same with the intpoints() option.

2. When using the RR after nlcom with both xtlogit and meglm, the p-values seem to be incorrect. . . . This p-value seems too small compared to the confidence interval. For one of the sub-group models the 95% CI overlaps 1, while the p-value is still very small: . . . Any ideas why this might be and how to get the correct p-values for RRs? The p-values for the RDs seem to be ok.

Discussed above in the thread in #7. Unless the null hypothesis that the quantity estimated by nlcom equals zero is what you're looking for, as in RD, then I would just ignore all of the default p-value stuff spewed out by nlcom just as for margins.

3. When testing for interaction between arm and other variables using LRT, these are calculated for the underlying model, not the RR derived from nlcom. For mixed-effects logistic models most interactions appear to be highly significant because of the strong effects sizes with ORs, whereas the effect sizes and stratum-specific effects are much smaller when presenting RRs. So to get the RRs directly from the same models that are used to test for interactions, I tried xtpoisson . . . These both work ok and the results are similar to each other. Is it appropriate to use poisson models with these data, and if so, which is better?

There's been literature over the past couple of decades or so, especially it seems in the epidemiology field, that advocates Poisson regression (with robust standard errors) to estimate RR in cohort studies and randomized controlled trials. I'm not especially familiar with it, and I'm not sure how much of what's been looked at extends cleanly to hierarchical / mixed-effects regression models. In what of this literature that has dealt specifically with cluster-randomized trials, the emphasis seems to have been on marginal models, i.e., GEE with cluster-robust standard errors.

More to your point: at the beginning, you mention

The preferred option from your suggestions was mixed-effects logistic with xtlogit, as this fits best with our statistical analysis plan.

As a principle, I'd try to stick with the plan and avoid post hoc changes to your study's protocol, especially, changes to the statistical method that are made on the basis of what's seen in the study's results.

4. I also tried the linear regression (regress) and random effects linear regression (xtreg) models with nlcom, but . . . don't run for some sub-group analyses.

These are both noniterative (least-squares) methods and should at least give results regardless of what you throw at them. By "don't run", are you saying that they omit some categories as collinear, have missing values for standard errors or give inflated, unrealistic standard errors?

Any advice on which model is most appropriate and which effect measure is most appropriate with these data would be much appreciated!

If it were me, I'd try to stay with what's prescribed in the statistical analysis plan, which seems to have called for fitting a hierarchical / mixed-effects logistic regression model. If your audience expects something different, e.g., RR, then as a first attempt, I would try to derive it from the prescribed model.

If you have specific concerns, for example, about whether non-collapsibility of ORs affects covariate-adjusted results (your third item above), then perhaps a more direct estimation of RR could be justified on the basis that these are exploratory secondary analyses. (I take it that your protocol didn't call for stratified randomization on these covariates, teenage category, district and whatnot.)

[Sonia Lewis]

Hi Joseph,

Thanks for the feedback and explanation of how to get the correct p-values from -nlcom-. This all seems to be fine now, but a colleague said that the 95%CIs of RD/RD using -nlcom- might be not correct because they may ignore the correlation between the estimated risks and the skewness in the distribution of RR/RD. Do you have any thoughts about this? I tried comparing confidence intervals for RR produced by models and produced by -nlcom- using these two approaches with poisson models:

*RR from poisson model
meglm abu i.armx i.age_cat16x districtx || locationx: , family(poisson) link(log) eform

*RR from -nlcom- after poisson model
quietly margins armx, post
nlcom (RR: _b[1.armx] / _b[0.armx]), noheader post
test RR = 1

I don't know if this exactly tests the concern my colleague raised, but I can see that the effect sizes are exactly the same, and there are indeed some small differences in 95% CIs between the two approaches. The differences are small enough that it does not affect interpretation. Is this anything to worry about? What about when the same is applied to the logistic models?

[Joseph Coveney]

. . . a colleague said that the 95%CIs of RD/RD using -nlcom- might be not correct because they may ignore the correlation between the estimated risks and the skewness in the distribution of RR/RD.

If there's a question as to the coverage of your risk-ratio or risk-difference confidence interval (or any other confidence interval), then you can determine its coverage via simulation.

In addition, it's not necessary to use -nlcom- to get the confidence interval for the risk-ratio or risk difference; rather, you can bootstrap the ratio or difference and use one or another of the bootstrap confidence intervals.

I tried comparing confidence intervals for RR produced by models and produced by -nlcom- using these two approaches with poisson models:[code elided for brevity] . . . I can see that the effect sizes are exactly the same, and there are indeed some small differences in 95% CIs between the two approaches. The differences are small enough that it does not affect interpretation. Is this anything to worry about? What about when the same is applied to the logistic models?

Again, I'm not sure whether the use of Poisson regression to estimate risk ratios extends cleanly to mulitlevel / hierarchical models (although the literature seems to have accepted the use of GEE for longitudinal data), but if the difference between the various approaches is small enough to be inconsequential (see the bottom-line summary in #2 above), then I wouldn't let it keep me awake at night. I would nevertheless strive to adhere to the study's protocol (statistical analysis section) as closely as feasible, absent any good reason to deviate from it.