nlcom gives very small p values for multiplicative interaction

[Mustafa Hussein]

Dear all,

I am estimating the association of socioeconomic status (low, middles, vs. high SES) with incident cardiovascular disease (cvd=1 if participant had an event, 0 otherwise) using a Poisson model with an offset for person-years in Stata 13.1. My sample size is > 5,000. I am interested in testing whether SES modifies the association of incidence with depressive symptoms (a continuous z score of the Centers for Epidemiologic Study Depression (CES-D) scale). Therefore, I have an interaction term of SES with zCES-D. This interaction term in the main model results provides a direct assessment of multiplicative interaction of the two variables. I also use nlcom to assess both additive and multiplicative interactions, following margins estimation. The results of nlcom and coefficient interactions from main model, although virtually identical in terms of point estimates, SEs, and CIs, are wildly different in the z statistics and p values. In models adjusted for further sets of covariates, nlcom starts using t statistics and the results remain quite divergent, in terms of statistical significance, from what the main model estimates show. Further, nlcom seems to almost always give very large z (or t) statistics (and very small p values) when testing multiplicative interactions, regardless of what the main model estimate and confidence intervals indicate.

Clearly, in this example (and others I have similarly tested), there is no interaction; the multiplicative-scale relative excess risk due to interaction (RERI) is too small and the confidence interval spans 1, and z statistics in the original model results support this. Why are nlcom z values (and p values) discrepant? I tried changing the degrees of freedom option of nlcom but it did not help (even unrealistically setting it to zero!).

Below is an example of my code and output. Compare estimates shown in boldface.

poisson cvd b3.ses##c.zcesd agecat2 fem2 race2 race3 race4 if sample==1, offset(log_cvdyr) irr

------------------------------------------------------------------------------
cvd | IRR Std. Err. z P>|z| [95% Conf. Interval]
-------------+----------------------------------------------------------------
ses |
Lower | 1.729769 .219026 4.33 0.000 1.349608 2.217015
Middle | 1.308588 .1600963 2.20 0.028 1.029591 1.663188
|
zcesd | .991996 .1206691 -0.07 0.947 .7815692 1.259077
|
ses#c.zcesd |
Lower | 1.089597 .1439638 0.65 0.516 .8410084 1.411665
Middle | 1.308522 .1835801 1.92 0.055 .9939425 1.722666
|
agecat2 | 2.640311 .2270739 11.29 0.000 2.230742 3.125077
fem2 | .4832043 .0411882 -8.53 0.000 .40886 .5710668
race2 | .5720179 .0882649 -3.62 0.000 .4227323 .774023
race3 | .7737112 .0816596 -2.43 0.015 .6291305 .9515181
race4 | .8091268 .0929864 -1.84 0.065 .645944 1.013534
_cons | .0070206 .0008026 -43.38 0.000 .0056113 .0087838
log_cvdyr | 1 (offset)
------------------------------------------------------------------------------

margins, at(zcesd=(0 1)) over(ses) expression(1000*(predict(ir))) post

1._at : 1.ses
zcesd = 0
2.ses
zcesd = 0
3.ses
zcesd = 0

2._at : 1.ses
zcesd = 1
2.ses
zcesd = 1
3.ses
zcesd = 1

------------------------------------------------------------------------------
| Delta-method
| Margin Std. Err. z P>|z| [95% Conf. Interval]
-------------+----------------------------------------------------------------
_at#ses |
1#Lower | 12.48939 .7401571 16.87 0.000 11.0387 13.94007
1#Middle | 9.78306 .699081 13.99 0.000 8.412887 11.15323
1#Upper | 7.625814 .7457161 10.23 0.000 6.164237 9.087391
2#Lower | 13.49948 .9569949 14.11 0.000 11.6238 15.37515
2#Middle | 12.69889 1.244115 10.21 0.000 10.26047 15.13731
2#Upper | 7.564777 1.393971 5.43 0.000 4.832644 10.29691
------------------------------------------------------------------------------

nlcom (lowhi_mRERI: (_b[2._at#1bn.ses]/_b[1bn._at#1bn.ses]) / (_b[2._at#3.ses]/_b[1bn._at#3.ses]) ) ///
(midhi_mRERI: (_b[2._at#2.ses]/_b[1bn._at#2.ses]) / (_b[2._at#3.ses]/_b[1bn._at#3.ses]) ), post

------------------------------------------------------------------------------
| Coef. Std. Err. z P>|z| [95% Conf. Interval]
-------------+----------------------------------------------------------------
lowhi_mRERI | 1.089597 .1439638 7.57 0.000 .8074333 1.371761
midhi_mRERI | 1.308522 .1835801 7.13 0.000 .9487121 1.668333
------------------------------------------------------------------------------

Thanks for your attention.

Mustafa

Mustafa Hussein, PhD
Postdoctoral Fellow
Drexel University School of Public Health
Philadelphia, PA, USA

[Joseph Coveney]

I thought that nlcom was supposed to be run after an estimation command, such as those that fit regression models and that have erturn matrices, scalars and the like. It seems that your having nlcom use coefficients returned by margins, which isn't one of those.

You might want to see whether it helps to have nlcom reference the regression coefficients of the fitted model returned by the poisson estimation command, itself, instead of those returned by margins.

[Joseph Coveney]

I take that back. It seems that the post makes margins act as an estimation command. So that can't be the reason.

[Maarten Buis]

nlcom uses an approximation (the delta method) to get at an estimate of the sampling distribution. This approximation often does not work well for ratios of coefficients, which is what you have found.

[Mustafa Hussein]

Thanks Maarten and Joseph for chiming in. Would not the use of the delta method also/rather affect the standard errors (SEs)? The SEs from the model estimates and from nlcom are identical, but interestingly not the confidence intervals, which are still close but not identical. I'd imagine once one has the SE and point estimate, the calculation of the z or t statistic should be straight forward?

[Joseph Coveney]

I believe that the Wald test statistic reported in the regression output table for poisson is computed from the untransformed coefficients and their asymptotic covariance matrix. You cannot really use the displayed IRR coefficient and its displayed standard error to compute the Wald test statistic shown next to it in the regression coefficient output table. Try it:

.ÿdisplayÿinÿsmclÿasÿtextÿ1.089597ÿ/.1439638
7.5685485

That matches what nlcom gets, but that's not what is shown in the regression output table (0.65), which is computed from the native regression coefficient and its standard error. The IRRs shown there (and their standard errors, which are also computed, I believe, with the delta method) are for display (as an aid in interpretation). I think, as Maarten mentions, you really cannot rely on the delta method for standard errors of ratios (and thus their test statistics) in these nonlinear models.

[Maarten Buis]

I would be careful with using nlcom for ratios, but the real reason for the difference in p-values is that they refer to different null hypotheses: with poisson the null hypothesis is that the ratio equals 1 (no effect) while with nlcom the null hypothesis is that the ratio equals 0 (usually not the hypothesis you are interested in for a ratio).

You can use the output of poisson to get close to the p-values: di (1.089597-1)/.1439638. So we create the z-value for the null-hypothesis that the incidence rate ratio equals 1. This produces output that is close but not quite the same as the z-value reported by poisson, because poisson uses the untransformed coefficients and standard errors.

[Mustafa Hussein]

Thank you guys so much for the input. Following nlcom, testing whether the calculated ratios are different from 1 gives close p values (based on Chi-sq distribution) to the Poisson model, which supports what Maarten mentioned about the tested null hypothesis.

Now, nlcom can be misleading if one does not pay attention to what it is doing. Is there any alternative to nlcom for calculating ratios and their inference statistics? Several interaction analysis tutorials recommend using nlcom without alerting people to the fact that nlcom is testing the null hypothesis the the ratio is equal to zero.

[Mustafa Hussein]

And, Maarten, I just found your Stata tip 116, which also touches on this issue with nlcom. Thanks for this decent tip.