Cox regression assessing joint relationship between baseline groups? JAMA example

Question

I'm a new R user trying to better understand the analytical framework behind an important JAMA paper (doi:10.1001/jama.291.2.210) and how it can be coded in R. It's objective was to evaluate if, in asymptomatic individuals, Coronary Arterial Calcium Score (CACS) assessment combined to the Framingham Risk Score (FRS) could provide prognostic information superior to either method alone. I'm using the survival package.

I'm not including a reproducible example as I understand this has to be figured out theoretically first.

1. Univariate Cox Model

For evaluating their prognostic value alone, two separate univariate Cox regressions were conducted. For the FRS model, they stratified individuals in four groups based on baseline FRS: 0-9%, 10-15%, 16-20% and 21% or higher. Let's call these 'FRS', with f1 f2 f3 f4 groups. For the CACS model, they stratified individuals based on baseline CACS: 0, 1-100, 101-300, 300 or higher. Let's call this CACS, with c1 c2 c3 c4 group.

This is straightforward, absolutely no problem where. Easily codable. HRs of 1.00 for reference hazards (c1 and f1).

2. Bivariate Cox Model - where I'm stuck

For evaluating joint effect of CACS and FRS over prognosis, they ran "bivariate cox regression models". What came out of this paper was a major contribution that changed clinical practice in primary prevention in Cardiology. Results below:

I understand this could elementary for some users, but I am doing a bit of homework and sounds like an interesting discussion.

Let's go:

First of all, I understand there are big terminology issues, as well addressed in here. A univariate Cox regression is a model considering only time-to-event and grouping variable. A Cox multiple regression is a model which include more than one predictor (often a given grouping factor and relevant covariates - often suboptimally called multivariate Cox regression). A true multivariate Cox regression would model time-to-event for more than one outcome with predictors. Bivariate cox regression is applied in cases where, due to impossibility to rule out one event as the outcome, such as recurring (e.g. acute manifestations of chronic diseases) or parallel events (e.g. retinopathy in left and right eye), event and trajectory to event have to be considered.

Cox and Oakes, "Analysis of Survival Data" 1982, Chapter 10, "Bivariate Survivor Functions":

This chapter is mainly concerned with applications which do not involve the singling out of one variable as a response, that is with studies of correlation rather than of regression. Then the full joint distribution of and becomes of interest. Consideration of the process unfolding in time directs attention to the hazard functions.

Q1: is the the model in the JAMA paper a bivariate cox analysis?

It doesn't seem so, as they don't focus on recurring myocardial infarctions. The outcome is non-fatal MI or CHD death. I understand they used BIvariate because of the joint relationship between predictors, but still, doesn't strike me as appropriate, since literature on Bivariate analysis usually analyze joint impact of time of recurring events. I don't hold any illusions that big journals only use right terms, but IDK, maybe one could argue this is a bivariate analysis.

Q2: what type of analysis is this?

It seems that this is a stratified Cox regression. Some posts cover more the stratification of covariates, not of the treatment/baseline group.

I'm too much of a beginner, but the logical way to go about it appears be to manually stratifying (could do within R or in a spreadsheet app) groups. No strata() coding. Just creating a new column (CACS_FRS) for defining groups combining FRS and CACS. In this case of CACS and FRS, it yields 16 groups. Then, selecting a reference for hazard function would be intuitive: c1f1 (CACS 0, FRS 0-9%). Would also include covariates of interest if I were to apply this model to other analyses. Does this idea seem appropriate? Would this be a valid way to evaluate joint relationship of grouping factors over prognosis?

xxxx$CACS_FRS <- factor(xxxx$CACS_FRS, levels = c(c1f1, ........))
coxph(Surv(time,event) ~ CACS_FRS + covariate1 + ... + covariaten, data = xxxx)

Q3: can you explain the existence of 4 reference HRs in the "Table 3 - Bivariate Analysis"?

If the way I mentioned above makes some sense, it's not what they did. I can't possibly come up with an explanation for 4 reference HR's. It would make some sense if they decided to stratify within each category, like: a model FRS 0-9% with the four CACS groups. But still, they would need either column 1 or line 1 to contain only HR = 1.00.

They state: The referent group has the lowest FRS (0%-9%) and a CACS of 300 or less or a low intermediate FRS (10%- 15%) and a CACS of zero. These groups were chosen as the referents due to simi- lar event rates.

Q4: how would you implement this type analysis?

I would really like to hear any tips or comments on how to do this, it just got me thinking and researching for some days now. Let me know if you guys have any R tips to code this - if it's more complex than it seems or what. Packages, functions, etc.

I found some other papers now that describe something like what I mentioned in my comment on Q2, like this one: doi:10.1001/jama.292.10.1188 .

Thank you very much!

Answer 1

This is not a bivariate Cox model. It is a Cox model with two predictors ("multiple" not "multivariate"). "Bivariate" is an unfortunate use of terminology.

Nor is it even a stratified Cox model in the sense of having more than one baseline hazard function-- "stratified" here means that you look at associations with one factor separately according to levels of another.

There are four reference HRs because the reference group is those four categories -- the paper say "These groups were chosen as the referents due to similar event rates". Every other individual cell in the table is compared to the combination of those four cells. I don't think that's good practice, but it does let you get bigger (but still finite) HRs elsewhere in the table.

How to implement? Take a complete cross-classification of Framingham Risk Score and CACS risk score to get a 16-level variable identifying the cells of the table. Now collapse the four reference levels together. For example, if the two scores were coded 1-4, then in R

cell<-paste(CACS,FRS,sep=":")
cell[cell %in% c("1:1","2:1","3:1","1:1")]<-"ref"
coxph(Surv(ttohardchd, hardchd)~factor(cell))

where hardchd is the indicator for coronary death or non-fatal MI, and ttohardchd is the corresponding observation time.

Update: So what model would make more sense? It's not feasible to use the upper left cell as the reference group, because there are no events in that group and the maximum partial likelihood estimator for the HRs in other cells will all be infinite. I might pick c2f2, or c2f3 (which is the median group on each margin) and then use the 16 cells as indicators. But I'd be more likely to try to model FRS and CACS as continuous variables.