For this example, we will use Approach B to conduct an interaction analysis based on Cox regression. We want to see if sex (z) moderates the association between body mass index (x) and out-patient care due to CVD (y).

Note Here, we use the same stset as in Cox regression.

Dataset

StataData1.dta

Variable name	cvd
Variable label	Out-patient care due to CVD (Ages 41-50, Year 2011-2020)
Value labels	0=No 1=Yes

Variable name	bmi
Variable label	N/A
Value labels	N/A

Variable name	sex
Variable label	Sex
Value labels	0=Man 1=Woman

Define the analytical sample

We start by defining the analytical sample:

gen pop_interact4=1 if cvd!=. & bmi!=. & sex!=.

Let us have a quick look at the variables:

sum cvd bmi sex if pop_interact4==1

Simple regression models

First, we will run the simple models, one for bmi and cvd, and one for sex and cvd.

stcox bmi if pop_interact4==1

stcox sex if pop_interact4==1

The tables above show that there is no association between bmi and cvd, whereas there is a clear and statistically significant association between sex and cvd.

Multiple regression model

Next, we run a model with both independent variables included:

stcox bmi sex if pop_interact4==1

The HRs for bmi and sex have decreased a tiny bit (in this case, become further from 1).  

Now, we need to save the estimates from this model:

estimates store model1

Multiple regression model with interaction effect

In this step, we will include the interaction term using Approach 2 (two hashtags mean that we specify the main effects and the interaction effect at the same time):

stcox c.bmi##i.sex if pop_interact4==1

Similar to the previous model, we have to save the estimates from this one:

estimates store model2

Note Since bmi is continuous and sex is binary, we actually do not need to compare the model fit using Approach B – it would be sufficient with Approach A (we can already see in the table that the interaction term is statistically significant). But since Approach B is what the example is about, we will do it for the sake of practice.

Next, we compare the fit of the two models:

lrtest model1 model2, stats

Note We called our saved models “model1” and “model2”, but you can choose any name you like.

In the table above, we can see that the p-value for the likelihood ratio test is below 0.05 (0.0225), which suggests that model that contains that interaction term (model2) fits the data better than the model without the interaction term (model1). This is also confirmed by the values for AIC, which are lower for model2. The BIC value is not lower for model2 than model1, but it should be noted that BIC tend to penalise complex models more than AIC does. 

We can thereby conclude that there is a statistically significant interaction between body mass index and sex in the effect on out-patient care due to CVD.

Interpretation

How can we understand this interaction effect that we found? Since our x-variable – sex – is binary, the easiest strategy for gaining more insight would be to do sex-specific analyses of the association between bmi and cvd. 

We start with a model for men, and then continue with the same for women.

stcox bmi if pop_interact4==1 & sex==0

stcox bmi if pop_interact4==1 & sex==1

The sex-specific models show that the association between body mass index and out-patient care due to CVD is actually negative (and statistically significant) for men and positive (and statistically non-significant) for women. This explains why we did not find any association at all in our first simple model.

Illustration

We would also like to show how the interaction effect can be illustrated. Unfortunately, it not very straight-forward to graph interactions for non-linear outcomes. A solution here might be is to categorise bmi. For example, we could categorise it into underweight, normal weight, overweight, and obese in according to WHO standards.

gen bmi_cat=bmi

recode bmi_cat (0/18.4999=1) (18.5000/24.999=2) (25.000/29.999=3) (30.000/40=4)

Let us also create and apply some value labels:

label define bmi_cat 1 "Underweight" 2 "Normal weight" 3 "Overweight" 4 "Obese"

label values bmi_cat bmi_cat

We choose Normal weight (value 2) as out reference category for bmi, and Women (value 1) as our reference category for sex.

Note Even though we usually do not have to specify a reference category for binary variables (such as sex), margins will not work without it.

First, we run a Cox regression model with these two variables. We use the quietly option to suppress the output.

quietly stcox ib2.bmi_cat##ib1.sex if pop_interact4==1

Then we type the margins command, with suppressed output here as well:

quietly margins ib2.bmi_cat#ib1.sex

Note We do not use double hashtags here to specify the main effects and the interaction effect, since it produces an additional line in the marginsplot.

And run the marginsplot:

marginsplot

Note The reference combination (predicted HR=1) includes women who are normal weight. All other estimates are compared to this combination. If we were to choose another reference combination, the graph would look quite different.

This graph illustrates the negative association for men and the positive association for women. Another way of looking at it is that the difference between men and women when it comes to the association between bmi_cat and cvd is largest among the underweight but becomes smaller as body mass index increases. Generally, the results might not be what we would expect (but remember that this dataset is fictional).

Summary Gender moderates the association between body mass index and out-patient care due to CVD. More specifically, a higher body mass index is associated with a lower risk of experiencing out-patient care due to CVD among men. Among women, the opposite association was observed.