# Is there a smart way to combine Cox survival analysis on a nested experimental design?

I have a study that ran a fairly simple survival experiment (control group, treatmentA and treatmentB), but was run under 3 different temperature conditions (10C, 30C, 35C). So i have 9 groups of individuals in total. I can run a survival analysis (coxph in R) within a particular temperature to see if my treatments influence survival, and i can run a survival analysis within a treatment type comparing across temperatures to see how temperature may influence survival.

What i would like to do is look overall at all 9 groups together (6 if you dont count the control/reference groups) rather than just 3 at once.

My limited experience brought me to thinking about running the 3 separate 'within temperature' models then pulling out the hazard ratios for each temp*treatment and comparing them visually - but are there better ways to do this?

When i look at the effect of temperature on the 3 control groups, it does seem to make a difference so i guess im needing to make sure each temperature control group is used as a reference for the corresponding treatment groups.

Any ideas of things to try or read would be welcome.

i considered using something like metafor to effectively treat it as a meta analysis of 3 studies but that seem to use the global outputs rather than the HRs from each treatment.

The trick here is to include interaction terms that cover each combination of treatment and temperature. You have to be careful in interpreting the coefficients returned by the analysis, but you will generally do much better with a combined analysis, taking advantage of all your data, than with separate analyses.

In R, you could simply write the model as:

coxph(Surv(time, status) ~ treatment * temperature


and the software will automatically generate all the implied coefficients. Specify the control group as the reference for the treatment categorical predictor. For temperature, code it as a 3-level categorical predictor, or keep it as numeric if you think that the log-hazard is linearly related to temperature.

If you code temperature as categorical with 10C as the reference, you will get the following coefficients:

treatmentA and treatment B are their log-hazard differences from control at 10C.

temperature30C and temperature35C are the log-hazard differences from 10C under control conditions.

The interaction coefficients treatmentA:temperature30C, treatmentA:temperature35C,treatmentB:temperature30C, and treatmentB:temperature35C represent the further log-hazard differences under those conditions from what you would get by simply adding the corresponding treatment and temperature coefficients.

So that in principle gives you all 8 differences in log-hazard (or hazard ratios, after exponentiation) from the baseline situation (control at 10C), by combining the coefficient values and their estimated errors appropriately. The emmeans and rms packages in R can help with converting the regression coefficients into estimates for specific conditions.

If you can treat temperature as a continuous numeric predictor linearly related to log-hazard things might be a bit simpler. To simplify interpreting the coefficients, subtract 10 from all of your temperatures, because by default the baseline conditions are reported with continuous predictors set at values of 0.

With a linear treatment of temperature, the coefficients for treatmentA and treatment B would be interpreted as above (at 10C). There would be a single temperature coefficient representing the change in log-hazard per degree under control conditions, and two interaction coefficients, treatmentA:temperature and treatmentB:temperature, representing the differences from the control temperature coefficient under those treatments.

Including interactions should work well if your study has enough events. You can evaluate the significance of the interaction terms; if they aren't significant then there's no evidence that the effect of treatment differs with temperature (or vice-versa) and you might simplify the model with simple additive treatment and temperature terms.

There's nothing specific to Cox regression modeling here except for the baseline hazard conceptually replacing the intercept of a linear regression model and the requirement to check that the proportional hazards assumption is met. So if you are confused by how to interpret coefficients when there are interaction terms, any reputable statistics text or many questions taged interaction on this site should help.