Cox Regression: Testing for effect in subgroup

Question

I'm using a cox proportional hazard model in R to see if a treatment variable (treatment or placebo) has effect on the survivaltime of patients. I intend to test this for each of my grouping variables (e.g. Age <60 or Age > 60).

I could do this by making two different subgroups of my dataset, as follows:

c.example1 <- coxph(Surv(time, status)~factor(treatment), data=lung[ lung$age.ind==0 ,])
c.example2 <- coxph(Surv(time, status)~factor(treatment), data=lung[ lung$age.ind==1 ,])

Each cox model gives me a P value telling me if there is a significant treatment effect. So far so good.

Now someone working with me on this problem suggested to use an interaction term to test the same thing. He feels this also gives answer to the same question, but requires less work (instead of two models, we build one):

c2.2 <- coxph(s ~ age + treatment + age:treatment, data=lung)

His idea is that the p value for our interaction term now tells us that if there is a significant effect for treatment in both subgroups. I thought it would only give information if the effect of age on treatment is different in the two subgroups. Not test is the effect of treatment is significant.

Is his interaction term approach the right one? If not, what is a good approach (besides coding two seperate models)?

Answer 1

Your initial approach, breaking down the age variable into subgroups for separate analyses, loses a lot of potentially valuable information. See this page, or many others under the binning tag on this site, for why that isn't a good idea.

The single model that your colleague suggested accomplishes two things beyond your separate subgroup analysis.

First, and most important for your original question, consider the model s ~ age + treatment, without the interaction term. A significant coefficient for treatment in this case would mean significance regardless of age. So that part of your colleague's suggestion does directly test whether the effect of treatment is significant. That would be the case whether you persist in grouping by age or use age as a continuous variable.

Second, in the model suggested by your colleague, the added interaction term further evaluates whether there is a difference in treatment effects as age changes, as your initial question and the answer from @DWin recognize. That may be very important to know, although you typically need more data to identify a significant interaction term than you need for a main effect.

So your colleague's suggested approach simply and efficiently answers your original question and might gather more useful information.

You will, of course, have to check that the proportional hazards model adequately describes your data, or else none of these models will provide the answers you seek.

Answer 2

As far as R programming goes, the formula you offer is equivalent to:

s ~ age * treatment 

The asterisk-operator is overloaded in R. When its arguments are numeric, it is multiplication, and when they are tokens or symbols in an R formula they are interactions that get expanded to all first order and second order ( and higher if there are more than two symbols.) Your representation of that model is also syntactically legal. We would have more to talk about if you had presented the output. I originally thought your were using the 'lung' dataset from pkg:survival, but it has no treatment variable.

It's true that SO tries to get interpretation and methodology questions shunted over to CrossValidated.com, but since that may happen anyway, I'll start that process. The interaction model gives you more information. Your question would be answered by a comparison of two models:

   s ~ age * treatment   # three coefficients
   s ~ age + treatment   # two effect coefficients

When there is no major difference in the "age effect" there will be no material change in the estimated coefficient for treatment across a comparison of those two models. The estimated "effect" in the "0"-groups will not be reflected in the treatment coefficients, since R uses treatment contrasts. The estimated effects in the older age group will be a combination of two or three coefficients depending on which treatment group is under consideration.

The question of whether there is a significant difference in the effect of treatment in the two age groups will be reflected in the magnitude of the age:treatment interaction and the associated p-value in the output from the first of those models. The results from the second of those models will return an treatment-estimate that is essentially a weighted average of the two models you were originally using. The mean of your two treatment estimates should be close to that estimate if the numbers of your cases in the age categories are similar.

You should consider looking at the possibility of modeling age as a continuous variable.