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

Question

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),
    family=binomial(link="cloglog"), nAGQ=7, data=grid3@data)

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.

Answer 1

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.