We discuss connections between the Cox proportional hazards model and Poisson generalized linear models as described in Whitehead (1980). We fit comparable models to a sample dataset using coxph()
, glm()
, phmm()
, and glmer()
and explore similarities.
We generate proportional hazards mixed model data.
library(phmm) n <- 50 # total sample size nclust <- 5 # number of clusters clusters <- rep(1:nclust,each=n/nclust) beta0 <- c(1,2) set.seed(13) Z <-cbind(Z1=sample(0:1,n,replace=TRUE), Z2=sample(0:1,n,replace=TRUE), Z3=sample(0:1,n,replace=TRUE)) b <- cbind(rep(rnorm(nclust), each=n/nclust), rep(rnorm(nclust), each=n/nclust)) Wb <- matrix(0,n,2) for( j in 1:2) Wb[,j] <- Z[,j]*b[,j] Wb <- apply(Wb,1,sum) T <- -log(runif(n,0,1))*exp(-Z[,c('Z1','Z2')]%*%beta0-Wb) C <- runif(n,0,1) time <- ifelse(T<C,T,C) event <- ifelse(T <= C,1,0) sum(event) phmmd <- data.frame(Z) phmmd$cluster <- clusters phmmd$time <- time phmmd$event <- event
fit.ph <- coxph(Surv(time, event) ~ Z1 + Z2, phmmd, method="breslow", x=TRUE, y=TRUE) summary(fit.ph)
# cox partial likelihood pl <- c() for(h in 1:length(eventtimes)){ js <- phmmd$time == eventtimes[h] # j star j <- phmmd$time >= eventtimes[h] if(sum(js) > 1) stop("tied event times") pl <- c(pl, exp(phmmd[js, "linear.predictors"])/ sum(exp(phmmd[j, "linear.predictors"]))) }
Next we create data to fit an auxilary Poisson model as described in Whitehead (1980) using the pseudoPoisPHMM()
function provided in the phmm
package. This function also extracts the linear predictors as estimated from the Cox PH model so that we can calculate likelihoods and degrees of freedom.
ppd <- as.data.frame(as.matrix(pseudoPoisPHMM(fit.ph))) # pois likelihood poisl <- c() eventtimes <- sort(phmmd$time[phmmd$event == 1]) for(h in 1:length(eventtimes)){ js <- ppd$time == eventtimes[h] & ppd$m >= 1 # j star j <- ppd$time == eventtimes[h] if(sum(js) > 1) stop("tied event times") poisl <- c(poisl, ppd[js, "N"]*exp(-1)*exp(ppd[js, "linear.predictors"])/ sum(ppd[j, "N"]*exp(ppd[j, "linear.predictors"]))) }
Poisson likelihood:
(coxph.pois.loglik = sum(log(poisl))) coxph.pois.loglik - fit.ph$loglik[2]
Poisson degrees of freedom
(coxph.pois.df = length(fit.ph$coef) + sum(phmmd$event))
We fit an auxiliary Poisson GLM and note that the parameter estimates for z1 and z2 are identical to the coxph() fit, and the likelihood and degrees of freedom are as expected.
ppd$t <- as.factor(ppd$time) fit.glm <- glm(m~-1+t+z1+z2+offset(log(N)), ppd, family=poisson) summary(fit.glm) cbind(coxph.coef = fit.ph$coef, glm.coef = coef(fit.glm)[c('z1', 'z2')]) cbind(coxph.pois.loglik, glm.loglik=logLik(fit.glm))
The additional parameter estimates correspond to the estimated log baseline hazard, which we verify using the basehaz()
function.
bh <- basehaz(fit.ph, centered = FALSE) cbind( coxph.bh.step = log(bh$hazard - c(0,bh$hazard[1:(length(bh$hazard)-1)]))[1:5], glm.bh.step = coef(fit.glm)[1:5] )
set.seed(20200316) fit.phmm <- phmm(Surv(time, event) ~ Z1 + Z2 + (Z1 + Z2|cluster), phmmd, Gbs = 100, Gbsvar = 1000, VARSTART = 1, NINIT = 10, MAXSTEP = 100, CONVERG=90) summary(fit.phmm)
ppd <- as.data.frame(as.matrix(pseudoPoisPHMM(fit.phmm))) poisl <- c() eventtimes <- sort(phmmd$time[phmmd$event == 1]) for(h in 1:length(eventtimes)){ js <- ppd$time == eventtimes[h] & ppd$m >= 1 # j star j <- ppd$time == eventtimes[h] if(sum(js) > 1) stop("tied event times") poisl <- c(poisl, ppd[js, "N"]*exp(-1)*exp(ppd[js, "linear.predictors"])/ sum(ppd[j, "N"]*exp(ppd[j, "linear.predictors"]))) }
Poisson likelihood:
phmm.pois.loglik = sum(log(poisl)) phmm.pois.loglik - fit.phmm$loglik[1]
Poisson degrees of freedom
# Poisson GLMM degrees of freedom length(unique(x$cluster)) * x$nrandom + x$nfixed traceHat(fit.phmm, "pseudoPois") # + 2*sum(phmmd$event)
We fit an auxiliary Poisson GLMM, although with a general variance-covariance matrix for the random effects (phmm() only fits models with diagonal variance-covariance matrix).
library(lme4) ppd$t <- as.factor(ppd$time) fit.lmer <- glmer(m~-1+t+z1+z2+ (z1+z2|cluster)+offset(log(N)), data=ppd, family=poisson, nAGQ=0) summary(fit.lmer)$coef fit.phmm$coef logLik(fit.lmer) phmm.pois.loglik - logLik(fit.lmer)[1]
cbind( phmm.bh.step = log(fit.phmm$lambda)[1:5], glm.bh.step = fixef(fit.lmer)[1:5] )
Any scripts or data that you put into this service are public.
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.