Cox PH model and regression with log-log link function: is there a connection?

In this post: interpreting estimates of cloglog logistic regression, I read that when interpreting the results of a regression with a log-log link function, the exponent of the estimated coefficient is equal to the hazard (= probability of mortality) per unit time. So to arrive at the probability of mortality per unit time, I take the exponent of the intercept and then the exponent of each coefficient*variable product, and add the results together?

On the other hand, to interpret a Cox proportional hazards model as a hazard per unit time, you multiply coefficient*variable for each coefficient, add everything up, exponentiate it, and multiply the result by the baseline baseline hazard function (if you can work out what the baseline hazard is).

These two models sound fairly different to me. However, I've been experimenting with them and I get suspiciously similar coefficients. Can anyone explain why?

Here is the code in R. Variables are continuous except for cop99, a factor with 7 levels; surv is a survival model built from survival time and status (1=event, 0=survived):

coxmodel <- cph(surv ~ cost_mean + elev_mean + popn_mean + cop99 + PAs_mean,
data = grid3@data, x = TRUE, y = TRUE, surv = TRUE, time.inc=1)
glmmodel <- glmer(status ~ cost_mean + elev_mean + popn_mean + PAs_mean + (1|cop99),


Here is the output for the Cox PH model:

                         Model Tests    Discrimination
Indexes
Obs     36918    LR chi2    1525.11    R2       0.056
Events   2222    d.f.            11    Dxy      0.462
Center 1.6031    Pr(> chi2)  0.0000    g        1.268
Score chi2 1268.99    gr       3.552
Pr(> chi2)  0.0000

Coef    S.E.   Wald Z Pr(>|Z|)
cost_mean -0.5778 0.0272 -21.21 <0.0001
elev_mean -0.5786 0.0345 -16.77 <0.0001
popn_mean  0.0729 0.0083   8.83 <0.0001
cop99=30   1.8168 0.3802   4.78 <0.0001
cop99=40   2.3828 0.4017   5.93 <0.0001
cop99=50   2.0343 0.3775   5.39 <0.0001
cop99=110  0.5465 0.4989   1.10 0.2733
cop99=160  0.0572 0.4474   0.13 0.8982
cop99=170 -0.2884 0.4563  -0.63 0.5273
cop99=190  2.0738 1.0607   1.96 0.0506
PAs_mean  -1.5656 0.1502 -10.42 <0.0001


and here is the summary of the glmer model:

Random effects:
Groups Name        Variance Std.Dev.
cop99  (Intercept) 0.9941   0.997
Number of obs: 36918, groups:  cop99, 8

Fixed effects:
Estimate Std. Error z value Pr(>|z|)
(Intercept) -3.84387    0.37838 -10.159  < 2e-16 ***
cost_mean   -0.56986    0.02713 -21.002  < 2e-16 ***
elev_mean   -0.57268    0.03444 -16.630  < 2e-16 ***
popn_mean    0.06933    0.00851   8.147 3.74e-16 ***
PAs_mean    -1.57468    0.15089 -10.436  < 2e-16 ***


..surely the coefficients for the 4 continuous variables are too similar for this to be a coincidence?

If the baseline hazard in the Cox model is 1, it'll have no effect on the hazard. Equally, exp(-3.84387) is insignificantly small.

One final puzzle: this all looks to me like the continuous variables should have a big effect, not just a 'significant' one. But the r-squared is about 0.06 in both models. Doesn't that mean that, on the contrary, the model doesn't explain much?

P.P.S. bearing in mind @EdM's advice in the answer below, I tried a binomial approach with R-INLA as follows:

lattice_temp <- poly2nb(grid3, row.names = grid3@data$$id2) nb2INLA(paste(getwd(), "/lattice.graph", sep = ""), lattice_temp) lattice.adj <- paste(getwd(), "/lattice.graph", sep = "") # INLA requires a dataframe and spatial lattice as input grid10 <- grid3@data survinla <- inla.surv(grid3@data$$def_mean, grid3@data\$status)
coxinla <- inla(survinla ~ 1 + cost_mean + elev_mean + popn_mean + cop99 +
f(id2, model = "bym", graph = lattice.adj, scale.model = TRUE),
data = grid10, family = "coxph",
control.compute=list(dic=TRUE, mlik=TRUE, waic=TRUE),
control.hazard = list(hyper = list(prec = list(param = c(0.001, 0.001)))))


...I ran this without the spatial term too. Adding the spatial term improves the marginal log-likelihood from -2167.46 to 890.66 and the WAIC score improves from 4376 to 3651. Maybe that's the way to go.

Having said that, some of the Cox PH model coefficients do have a low standard error so I wonder if it would still have some predictive value? Not that I know how to extrapolate to future years - which I should read up on (this is really a predictive and not an explanatory exercise) - but I do have yearly data for 2 of the covariates, so the prospect of a Cox model with time-dependent covariates is still a temptation...unless I can work them into a more standard model.

• I suspect that the low standard errors in the Cox PH model come from the (probably unwarranted) assumption of independence among events. Your point estimates of survival might still be OK, but their standard errors will be too low. A potential problem with predictions involving time-varying covariates is survivorship bias. You'll have to judge the risk of that for yourself, based on your understanding of the subject matter.
– EdM
Commented Jun 29, 2021 at 13:53

There is, as you suspect, a fundamental connection between a Cox model and a binomial generalized linear model with a complementary log-log (cloglog) link.* That is described exactly for interval-censored data, as explained on this page with a link to the literature for further reading. The binomial regression in that circumstance with cloglog link is actually called a "grouped proportional hazards model." The intercepts in that case are related to the differences in cumulative baseline hazard between the endpoints of the time intervals.

The above is for a binomial model that explicitly includes (discrete) times in the model. With data organized so that censored cases are omitted beyond the last observation time, a discrete-time binomial model handles censoring like a continuous-time survival model: cases only contribute to calculations at times when they are in the risk set.

It's not clear that your binomial model handled times to events. Perhaps there was a fixed study duration with all participants starting at the same time, so that each individual either had an event during the study or was censored at a survival time equal to the full study duration. That would be OK: you then have interval censoring with a single interval equal to study duration, and the explanation above for the cloglog link would hold. In general, though, if you don't model time properly in a binomial model you risk having problems if there's censoring.

Finally, don't perseverate over the apparently small pseudo-$$R^2$$ value of the Cox model. As Frank Harrell explains:

The only thing going against pseudo $$R^2$$ is the difficulty in interpreting its absolute value. But it is excellent for comparing two or more models, even though examining increases in LR [likelihood ratio] is better.

With respect to your survival model, the $$D_{xy}$$ value shows how well you can distinguish between cases. You can calculate the concordance $$C$$ from that, $$C= 0.5 + (D_{xy}/2)$$ or about 0.73 in your Cox cph model. That means for any pair of cases whose survival times can be compared, the model puts them in the correct order 73% of the time. As a next step, use the calibrate() function for cph objects to get an estimate of how well predicted and actual survival probabilities align at any specified survival time.

*Technically you shouldn't call that a "logistic regression," as that term should be reserved for a binomial regression with a logit link.

• Many thanks @EdM. I got into trouble with the function 'calibrate' and asked a StackOverflow question: link ...the subject of the model is a bit unusual. I have 2km grids, and the event 'status = 1' is where deforestation in 30m pixels (over 2000-20) reaches 10% of the maximum seen in my dataset. The 'time to event' is the mean year of deforestation. For all grid squares with status = 0, I set the time to event as 21 (but maybe that wasn't necessary?) Commented Jun 28, 2021 at 8:52
• @MichaelSmith you might need to reconsider your approach. Your "events" are unlikely to be independent among neighboring grid squares (as your current Cox and binomial models implicitly assume). There are ways to take geospatial correlations into account in regression modeling, but I don't have any experience with them. Take care of that problem before you worry much more about improving the current, independence-assuming, models.
– EdM
Commented Jun 28, 2021 at 13:10
• This is definitely a concern. I'm looking at models with a lattice-based spatial element using R-INLA (hence the interest in log-log models) to see if the AIC value is better with/without a spatial element. Some of the variables (cost-distance as a measure of accessibility, and a categorical variable of vegetation type) will have information on space but I'm not sure if that deals with the problem. I was hoping to do a cross-validation exercise ...perhaps that's easier without using a coxph approach? It's really a predictive (not explanatory) project. Commented Jun 28, 2021 at 13:42