[Joseph Coveney]

We can't report the the risk ratio directly from the model results because of the interaction.

I agree that risk difference is probably preferable, but before closing:

.ÿ
.ÿversionÿ15.1

.ÿ
.ÿclearÿ*

.ÿ
.ÿsetÿseedÿ`=strreverse("1496832")'

.ÿ
.ÿlocalÿn0ÿ=ÿceil(0.107299ÿ*ÿ(1ÿ-ÿ0.107299)ÿ/ÿ0.0153276ÿ/ÿ0.0153276)

.ÿlocalÿn1ÿ=ÿceil(0.022342ÿ*ÿ(1ÿ-ÿ0.022342)ÿ/ÿ0.0123942ÿ/ÿ0.0123942)

.ÿ
.ÿquietlyÿsetÿobsÿ3

.ÿgenerateÿbyteÿsevÿ=ÿ_n

.ÿlabelÿdefineÿSeverityÿ1ÿMildÿ2ÿModerateÿ3ÿSevere

.ÿlabelÿvaluesÿsevÿSeverity

.ÿ
.ÿquietlyÿexpandÿ2

.ÿbysortÿsev:ÿgenerateÿbyteÿawaÿ=ÿ_n

.ÿlabelÿdefineÿAwarenessÿ1ÿUnawareÿ2ÿAware

.ÿlabelÿvaluesÿawaÿAwareness

.ÿ
.ÿgenerateÿintÿcount1ÿ=ÿrbinomial(cond(awaÿ==ÿ1,ÿ`n0',ÿ`n1'),ÿ0.05)

.ÿquietlyÿreplaceÿcount1ÿ=ÿcond(awaÿ==ÿ1,ÿrbinomial(`n0',ÿ0.107299),ÿrbinomial(`n1',ÿ0.022342))ÿinÿ5/6

.ÿgenerateÿintÿcount0ÿ=ÿcond(awaÿ==ÿ1,ÿ`n0',ÿ`n1')ÿ-ÿcount1

.ÿ
.ÿquietlyÿreshapeÿlongÿcount,ÿi(sevÿawa)ÿj(use)

.ÿ
.ÿ*
.ÿ*ÿBeginÿhere
.ÿ*
.ÿquietlyÿglmÿuseÿi.sev##ib2.awaÿ[fweight=count],ÿfamily(binomial)ÿlink(log)ÿnolog

.ÿcontrastÿr.awa@sev,ÿeformÿnowald

Contrastsÿofÿmarginalÿlinearÿpredictions

Marginsÿÿÿÿÿÿ:ÿasbalanced

------------------------------------------------------------------------------
ÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿ|ÿÿÿÿÿexp(b)ÿÿÿStd.ÿErr.ÿÿÿÿÿ[95%ÿConf.ÿInterval]
-----------------------------+------------------------------------------------
ÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿawa@sevÿ|
ÿÿÿÿ(UnawareÿvsÿAware)ÿMildÿÿ|ÿÿÿ1.343546ÿÿÿ.6019039ÿÿÿÿÿÿ.5583636ÿÿÿÿ3.232867
(UnawareÿvsÿAware)ÿModerateÿÿ|ÿÿÿ.4506303ÿÿÿ.1544947ÿÿÿÿÿÿ.2301406ÿÿÿÿ.8823632
ÿÿ(UnawareÿvsÿAware)ÿSevereÿÿ|ÿÿÿ6.834559ÿÿÿ4.910415ÿÿÿÿÿÿ1.671664ÿÿÿÿ27.94293
------------------------------------------------------------------------------

.ÿ
.ÿ//ÿAndÿthen,ÿforÿriskÿdifferences
.ÿmarginsÿsev,ÿdydx(awa)

ConditionalÿmarginalÿeffectsÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿNumberÿofÿobsÿÿÿÿÿ=ÿÿÿÿÿÿ1,653
ModelÿVCEÿÿÿÿ:ÿOIM

Expressionÿÿÿ:ÿPredictedÿmeanÿuse,ÿpredict()
dy/dxÿw.r.t.ÿ:ÿ1.awa

------------------------------------------------------------------------------
ÿÿÿÿÿÿÿÿÿÿÿÿÿ|ÿÿÿÿÿÿÿÿÿÿÿÿDelta-method
ÿÿÿÿÿÿÿÿÿÿÿÿÿ|ÿÿÿÿÿÿdy/dxÿÿÿStd.ÿErr.ÿÿÿÿÿÿzÿÿÿÿP>|z|ÿÿÿÿÿ[95%ÿConf.ÿInterval]
-------------+----------------------------------------------------------------
1.awaÿÿÿÿÿÿÿÿ|
ÿÿÿÿÿÿÿÿÿsevÿ|
ÿÿÿÿÿÿÿMildÿÿ|ÿÿÿ.0144145ÿÿÿÿ.020285ÿÿÿÿÿ0.71ÿÿÿ0.477ÿÿÿÿ-.0253434ÿÿÿÿ.0541724
ÿÿÿModerateÿÿ|ÿÿ-.0537845ÿÿÿ.0268508ÿÿÿÿ-2.00ÿÿÿ0.045ÿÿÿÿÿ-.106411ÿÿÿ-.0011579
ÿÿÿÿÿSevereÿÿ|ÿÿÿ.0816022ÿÿÿ.0175592ÿÿÿÿÿ4.65ÿÿÿ0.000ÿÿÿÿÿ.0471868ÿÿÿÿ.1160177
-------------+----------------------------------------------------------------
2.awaÿÿÿÿÿÿÿÿ|ÿÿ(baseÿoutcome)
------------------------------------------------------------------------------
Note:ÿdy/dxÿforÿfactorÿlevelsÿisÿtheÿdiscreteÿchangeÿfromÿtheÿbaseÿlevel.

.ÿ
.ÿexit

endÿofÿdo-file


.

contrast also works with survey data.

[Molly Jeffery]

Joseph Coveney: thank you for the suggestion. GLM with the non-canonical link does not converge with the actual data. I'm unclear, though, why you think the margins results would be different after GLM vs. logistic regression (contrast is part of the margins commands). Does -contrast- not use delta method standard errors?



Clyde Schechter: Yes, I agree. It is important to present the confidence intervals and p-values rather than just the result of the hypothesis test against the arbitrary p-value of 0.05. That's why I wanted to find out why two algebraically equivalent tests gave very different answers. In the past when we've done this, the differences in the p-values have been minor between the two methods. In this case, a p-value of 0.0002 vs. 0.1845 is clearly demonstrating that there is a problem somewhere. I appreciate your help determining that the problem was coming from the delta-method standard errors producing negative values in the confidence intervals for the rarest outcomes. Although we won't use the bootstrapping method in the paper, the exercise of making it work was very useful for me.

[Joseph Coveney]

Does -contrast- not use delta method standard errors?

There is no transformation involved. And you won't have the problem with the confidence bounds lying outside the parameter space.

I wonder whether the fact that you couldn't attain convergence with the risk ratio model and the fact that the p-values are so discrepant aren't unrelated.

[Joseph Coveney]

I wonder whether the fact that you couldn't attain convergence with the risk ratio model and the fact that the p-values are so discrepant aren't unrelated.

You can get such discrepant p-values with nlcom even with a dataset that allows convergence under all three model paramaterizations.

version 15.1

clear *

set seed `=strreverse("20190507")'
quietly set obs 6

generate byte sev = mod(_n, 3) + 1
label define Severity 1 Mild 2 Moderate 3 Severe
label values sev Severity

bysort sev: generate byte awa = mod(_n, 2) + 1
label define Awareness 1 Unaware 2 Aware
label values awa Awareness

generate double pr = 0.1 + (sev < 3) * 0.02 - (awa == 2) * (0.08 - (sev < 3) * 0.06)

generate byte count1 = rbinomial(100, pr)
generate byte count0 = 100 - count1
quietly reshape long count, i(sev awa) j(use)
table sev awa [fweight=count], contents(mean use min pr) format(%04.2f)

quietly glm use i.sev##ib2.awa [fweight=count], family(binomial) link(identity) nolog
assert e(converged)
contrast r.awa@sev

margins sev#awa, post
testnl _b[3.sev#2.awa] / _b[3.sev#1.awa] = 1
nlcom _b[3.sev#1.awa] / _b[3.sev#2.awa]

quietly glm use i.sev##ib2.awa [fweight=count], family(binomial) link(log) nolog
assert e(converged)
contrast r.awa@sev, eform

quietly glm use i.sev##ib2.awa [fweight=count], family(binomial) link(logit) nolog
assert e(converged)
contrast r.awa@sev, eform

margins sev#awa, post
test 3.sev#2.awa = 3.sev#1.awa
testnl _b[3.sev#2.awa] / _b[3.sev#1.awa] = 1

lincom _b[3.sev#1.awa] - _b[3.sev#2.awa]
nlcom _b[3.sev#1.awa] / _b[3.sev#2.awa]

exit

The helpfiles for testnl and nlcom both emphasize that the p-values based upon the delta method are "an approximation appropriate in large samples." I'm still not sure why testnl is so far off in the survey's data (it's able to distinguish a similar ratio from 1 in the code above), but I guess that "large" for these circumstances means large.

[Molly Jeffery]

Interestingly, on the actual data, svy: poisson converges, gives almost identical margins estimates as the logistic regression, including the negative values for the lower CI limit on the rarest outcomes, and happily calculates the contrast ratio with positive CI limits (4.95 [1.48 to 16.64])

Honestly, I wouldn't call the sample here large. I'm working with a subpopulation of 1300, and some of the strata (stage x awareness) are quite small (around 100). I usually work with samples of millions of people, where I sometimes don't bother to calculate standard errors because everything has p<0.001. Yet reviewers still wanted corrections for multiple comparisons . . .

Thanks again for your interest.