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.