How to calculate P for interaction of subgroup analysis using STATA?

[Wei Liu]

Hello, I am a graduate, and I am woking on a project on nutritional epidemiology. Now, I need to calculate P value for interaction of subgruop analysis. For exmaple, the exposure variable is q5_adj_food (categorical varibale, 5 categories), and the outcome variable is cvd_mortality. Now, I have obtained the hazard ratios by sex with the following command (multivariable Cox proportional hazard regression model):

by sex, sort: stcox age_dhq sex ib(first).race4 ib(first).center ib(first).marital ib(first).educat asp hyperten_f diabetes_f ib(first).cig_stat ib(first).bmi_curc dt_alc_dhq physical_activity dt_kcal_dhq ib(first).q5_adj_food,nolog 

-> sex = 1
failure _d: cvd_mortality == 1
analysis time _t: mortality_py
note: sex omitted because of collinearity
Cox regression -- Breslow method for ties
No. of subjects = 41691	Number of obs = 41691
No. of failures = 3169
Time at risk = 559931.6603
	LR chi2(34) = 1850.30
Log likelihood = -29712.597	Prob > chi2 = 0.0000
_t Haz. Ratio Std. Err.	z	P>z [95% Conf. Interval]
age_dhq 1.116506 .0037055	33.21	0.000 1.109267 1.123792
sex 1 (omitted)
race4
2 1.144675 .1093609	1.41	0.157 .949204 1.3804
3 .9179986 .1394609	-0.56	0.573 .6815987 1.236389
4 .6967727 .1354061	-1.86	0.063 .4760738 1.019783
center
2 1.110001 .101349	1.14	0.253 .9281203 1.327525
3 1.485132 .306829	1.91	0.056 .9906183 2.226505
4 1.221528 .0953126	2.56	0.010 1.048302 1.423379
5 .7401891 .0533708	-4.17	0.000 .6426396 .8525462
6 .9422014 .0826673	-0.68	0.497 .7933424 1.118992
8 1.174933 .0917145	2.07	0.039 1.008252 1.36917
9 1.004963 .0868551	0.06	0.954 .8483671 1.190463
10 1.002671 .0774231	0.03	0.972 .8618494 1.166502
11 1.169032 .1456933	1.25	0.210 .9156792 1.492483
marital
2 1.358208 .1025161	4.06	0.000 1.171435 1.574759
3 1.382282 .0890831	5.02	0.000 1.21826 1.568388
4 1.158425 .214747	0.79	0.428 .8055147 1.665951
5 1.567762 .1354457	5.20	0.000 1.323553 1.857031
educat
6 .8553797 .0419344	-3.19	0.001 .777015 .9416478
7 .8309011 .0406846	-3.78	0.000 .7548675 .9145931
asp .9743588 .0351426	-0.72	0.471 .9078587 1.04573
hyperten_f 1.533344 .0572585	11.45	0.000 1.425128 1.649778
diabetes_f 1.472549 .0823065	6.92	0.000 1.319754 1.643035
cig_stat
1 1.714912 .1009937	9.16	0.000 1.527965 1.924732
2 1.083727 .0434533	2.01	0.045 1.001821 1.17233
bmi_curc
2 .9716498 .3457569	-0.08	0.936 .483742 1.951667
3 1.095438 .3889924	0.26	0.797 .5461657 2.197108
4 1.371851 .4887097	0.89	0.375 .682455 2.757655
dt_alc_dhq 1.001684 .0005189	3.25	0.001 1.000668 1.002702
physical_activity .9994882 .0001497	-3.42	0.001 .9991949 .9997816
dt_kcal_dhq 1.000002 .0000276	0.08	0.937 .9999481 1.000056
q5_adj_food
2 .9434859 .0596177	-0.92	0.357 .8335833 1.067878
3 .9072255 .057343	-1.54	0.123 .8015181 1.026874
4 .9055417 .0591288	-1.52	0.129 .7967606 1.029175
5 .9650946 .0668769	-0.51	0.608 .8425298 1.105489
-> sex = 2
failure _d: cvd_mortality == 1
analysis time _t: mortality_py
note: sex omitted because of collinearity
Cox regression -- Breslow method for ties
No. of subjects = 48690	Number of obs = 48690
No. of failures = 2193
Time at risk = 655957.1562
	LR chi2(34) = 1943.98
Log likelihood = -20467.224	Prob > chi2 = 0.0000
_t Haz. Ratio Std. Err.	z	P>z [95% Conf. Interval]
age_dhq 1.152671 .0047998	34.12	0.000 1.143302 1.162117
sex 1 (omitted)
race4
2 1.058327 .1060052	0.57	0.571 .8696834 1.28789
3 1.028641 .2229481	0.13	0.896 .6726261 1.57309
4 .6221799 .1553546	-1.90	0.057 .3813959 1.014976
center
2 .9611455 .1222795	-0.31	0.755 .7490261 1.233336
3 1.325036 .3612462	1.03	0.302 .7765399 2.260955
4 1.280477 .1260641	2.51	0.012 1.055772 1.553006
5 .7877282 .0810268	-2.32	0.020 .6439033 .9636784
6 .788818 .0865119	-2.16	0.031 .6362432 .9779812
8 1.077666 .1123804	0.72	0.473 .8784554 1.322052
9 1.149042 .1181497	1.35	0.177 .9393155 1.405595
10 .8236091 .0905715	-1.76	0.078 .6639194 1.021708
11 1.115448 .1506414	0.81	0.419 .8560409 1.453464
marital
2 1.16392 .0645165	2.74	0.006 1.044096 1.297494
3 1.203423 .0836954	2.66	0.008 1.050072 1.379169
4 .9959414 .2689946	-0.02	0.988 .5865865 1.690968
5 1.292524 .1453415	2.28	0.022 1.036866 1.611218
educat
6 1.07853 .0675791	1.21	0.228 .9538874 1.219459
7 .8861971 .0603997	-1.77	0.076 .7753821 1.012849
asp .9980623 .0435925	-0.04	0.965 .9161775 1.087266
hyperten_f 1.591741 .0716107	10.33	0.000 1.457397 1.73847
diabetes_f 1.948909 .1338385	9.72	0.000 1.703478 2.229701
cig_stat
1 2.252348 .1543822	11.85	0.000 1.969209 2.576198
2 1.220426 .0585455	4.15	0.000 1.110909 1.340741
bmi_curc
2 .9300074 .1733491	-0.39	0.697 .6453961 1.340129
3 .8808197 .1651636	-0.68	0.499 .6099261 1.272029
4 1.028685 .195279	0.15	0.882 .7090805 1.492344
dt_alc_dhq .9995631 .0019559	-0.22	0.823 .9957368 1.003404
physical_activity .9992718 .0002112	-3.45	0.001 .9988579 .9996858
dt_kcal_dhq .9998937 .0000442	-2.41	0.016 .9998071 .9999802
q5_adj_food
2 .9889278 .0626882	-0.18	0.861 .8733872 1.119753
3 1.19714 .0796657	2.70	0.007 1.050752 1.363922
4 1.14269 .0828671	1.84	0.066 .9912878 1.317216
5 1.342852 .1025518	3.86	0.000 1.156174 1.559672

As you can see, the HRs are obviously different between men and women. Now, I want to use likelihood ratio test to test whether there is a statistical difference in HRs between two subgroups? In other words, I want to obtain the P value for interaction between my exposure variable and sex. Here, can you tell me the specific STATA procedure？Of note, my exposure of interest is a categorical varibale (5 categories), and the regression model is Cox proportional hazrd regression model.
Thanks for your time and help in advance.
Kind regards,
Dr. Liu

[Marcos Almeida]

Welcome to the Stata Forum / Statalist,

Please read the FAQ. There you'll learn how to share command/output/data in the Forum. Thanks.

That said, I gather you could: first, perform the model with sex (not by sex), hence having coefficients for the difference between sexes; applying interactions terms. For this you just need to add # between variables. If you want the so-called "main effects" plus interactions terms, the abbreviate command inserts ##. For example: stcox c.age_sq##i.sex

[Wei Liu]

Welcome to the Stata Forum / Statalist,

Please read the FAQ. There you'll learn how to share command/output/data in the Forum. Thanks.

That said, I gather you could: first, perform the model with sex (not by sex), hence having coefficients for the difference between sexes; applying interactions terms. For this you just need to add # between variables. If you want the so-called "main effects" plus interactions terms, the abbreviate command inserts ##. For example: stcox c.age_sq##i.sex

Dr. Marcos,
Thanks for your kindly reminder. I will read the FAQ latter. For your answer "perform the model with sex (not by sex)", do you mean that the Cox model only includes the variable sex (i.e., stcox i.sex, nolog)? If so, I feel it seems to be incorrect, as I want to control for other covariates. Can you state more specifically and clealy? Thanks you very much.
Best regards,
Wei Liu

[Marcos Almeida]

I meant: add sex as a predictor, plus the remaining predictors, plus modeling by adding (or not) interaction terms.

[Wei Liu]

I meant: add sex as a predictor, plus the remaining predictors, plus modeling by adding (or not) interaction terms.

Many thanks! I got it!

[Wei Liu]

I meant: add sex as a predictor, plus the remaining predictors, plus modeling by adding (or not) interaction terms.

Thanks again! Dr. Marcos, I have done it as you suggested. I am using stcox in Stata 12.0. The specific code is listed as below:

stcox age_dhq sex ib(first).race4 ib(first).center ib(first).marital ib(first).educat asp hyperten_f diabetes_f ib(first).cig_stat ib(first).bmi_curc dt_alc_dhq physical_activity dt_kcal_dhq i.q5_adj_food i.sex##i.q5_adj_food

.

The results are posted as below:

_t	Haz. Ratio	Std. Err.	z	P>z	[95% Conf.	Interval]
age_dhq	1.129674	.0029178	47.21	0.000	1.123969	1.135407
sex	.5298233	.0333472	-10.09	0.000	.4683346	.5993851
race4
2	1.095476	.0755315	1.32	0.186	.9570037	1.253984
3	.9441204	.1173633	-0.46	0.644	.7399712	1.204592
4	.677993	.1038063	-2.54	0.011	.502226	.9152741
center
2	1.058623	.0783953	0.77	0.442	.9156019	1.223986
3	1.414436	.2326298	2.11	0.035	1.024679	1.952446
4	1.261797	.0762124	3.85	0.000	1.120927	1.420372
5	.7530837	.0444192	-4.81	0.000	.6708675	.8453756
6	.877952	.0597071	-1.91	0.056	.7683921	1.003133
8	1.141434	.0711161	2.12	0.034	1.010224	1.289687
9	1.070988	.0694147	1.06	0.290	.9432243	1.216057
10	.9381442	.0591984	-1.01	0.312	.8290056	1.061651
11	1.139793	.1028131	1.45	0.147	.9550912	1.360214
marital
2	1.269783	.0562842	5.39	0.000	1.164124	1.385032
3	1.306919	.0617841	5.66	0.000	1.191265	1.433801
4	1.122243	.1713223	0.76	0.450	.8320361	1.513673
5	1.479129	.1010474	5.73	0.000	1.293766	1.691049
educat
6	.925967	.0357103	-1.99	0.046	.8585558	.9986711
7	.8463226	.0333812	-4.23	0.000	.7833617	.9143438
asp	.9816605	.0272823	-0.67	0.505	.9296185	1.036616
hyperten_f	1.559438	.0447541	15.48	0.000	1.474143	1.649669
diabetes_f	1.626498	.0703608	11.24	0.000	1.494278	1.770418
cig_stat
1	1.917528	.0854204	14.61	0.000	1.757207	2.092475
2	1.139902	.0348355	4.28	0.000	1.07363	1.210265
bmi_curc
2	.8817192	.1450609	-0.77	0.444	.63869	1.217224
3	.9215158	.1515971	-0.50	0.619	.6675334	1.272133
4	1.113961	.1845105	0.65	0.515	.8051591	1.541198
dt_alc_dhq	1.00187	.0004933	3.80	0.000	1.000904	1.002838
physical_activity	.9994263	.0001221	-4.70	0.000	.9991871	.9996656
dt_kcal_dhq	.9999767	.0000234	-1.00	0.319	.9999307	1.000023
q5_adj_ultrafood
2	.9564849	.0601905	-0.71	0.480	.8454988	1.08204
3	.9367556	.0584106	-1.05	0.295	.828992	1.058528
4	.9468759	.0603144	-0.86	0.391	.8357434	1.072786
5	1.040924	.0687684	0.61	0.544	.9145023	1.184823
2.sex	1	(omitted)
sex#q5_adj_food
2 2	.9931504	.0872765	-0.08	0.938	.8360126	1.179824
2 3	1.183747	.1030758	1.94	0.053	.9980207	1.404035
2 4	1.100383	.0992055	1.06	0.289	.9221543	1.313059
2 5	1.141487	.1035711	1.46	0.145	.9555169	1.363651

As you can see, Stata gave four P values (the red part). Now, I want to know how can I get an overall P for interaction?

Many thanks for your patience. I am looking forwar your reply.

Best regards,

Liu

[Marcos Almeida]

Just a clarifying note: if you use ##, you already included the so-called "main effects" plus the interaction term. There is no need to repeat it.

That said, if we have significant interaction, the best approach is providing a graphical display of margins, and that we get with - marginsplot - command. Please type - help margins - and help marginsplot - for that matter.

[Tala Smith]

Hello Wei,

I was wondering if you figured out how to determine the overall p-value for interaction and could share how you obtained it?

[Paul Dickman]

There's a section called "Testing significance of the interaction" on the page https://www.pauldickman.com/software...-interactions/ that illustrated two ways to determine the overall p-value for interaction.

[Tala Smith]

Would the methodology from the link you provided be the same for determining the p-value for interaction for subgroup analysis in a meta-analysis?