cloglog vs logistic regression survival analysis

Question

I'm not looking for an explanation of the difference between hazards and odds here, what I am curious about is determining if, as is sometimes presented, a hazard model (cloglog link) is equivalent to a binomial (logit link) when dealing with irregular interval censored data.

This comes up a lot in ecology for survival analyses, and for me in particular, nest survival. Not all nests are observed from their initiation, but ages are known, and visits are irregular, leading to interval censoring.

The more common approach is to use the logistic exposure method (power logistic; https://rpubs.com/bbolker/logregexp), but folks often perceive the hazard form (cloglog link) to be preferrable.

I have code to simulate and analyze a data set where survival increases with nest age, but I get different daily survival probabilities and cumulative survival probabilities between the two. I'm curious if I am dealing with the irregular interval censoring incorrectly somehow, back transforming predicted values wrong, or something else like the shape of the cdf of the cloglog vs logit distributions?

#simulate data
library(tidyverse)
set.seed(24)

max.age=25 #nest period length
N=40  #number of nests
age.seq<-seq(0,23,1) # age covariate values (easier to fix dsr @ age=1 at mu.s this way)
mu.s=qlogis(0.93) # logit survival probability, age=1
beta.age=0.3 #coefficent expressing the effect of nest age on dsr
S<-plogis(mu.s+beta.age*age.seq) # vector of dsr
cumprod(S) #cumulative survival probability


true.mat<-matrix(NA,nrow=N,ncol=max.age)

 true.mat[,1]<-1
 
 for(i in 1:N){
   for(t in 2:max.age){
     
     true.mat[i,t]<-true.mat[i,t-1]*rbinom(1,1,S[t-1])
   }
 }
 

true.mat%>%
  as.data.frame()%>%
  rownames_to_column(var="nest.id")%>%
  pivot_longer(cols=starts_with('V'),
               names_to = 'visit',
               values_to = 'ld')%>%
  mutate(is.obs=ifelse(visit=='V25',
                       1,
                       rbinom(n=nrow(.),size=1,prob=0.6) # probability a nest is visited on a given day
                       )
         )%>%
  filter(!(is.obs==0))%>%
  mutate( obs.dat=ld*is.obs)->nest.dat

#format data for analysis in glm, using just intervals where a nest was visited
# errors come up for nests that survived the entire period

nest.dat%>%
  mutate(day=as.numeric(gsub('V','',visit)))%>%
  select(-visit)%>%
  group_by(nest.id)%>%
  mutate(exposure=day-lag(day))%>%
  filter(!is.na(exposure))%>%
  filter(row_number()<=min(which(obs.dat==0))|
           obs.dat==1)%>%
  select(nest.id,day,exposure,
         obs.ld=obs.dat)%>%
ungroup()->glm.dat

Now for the analysis using the logistic exposure method. Code for the link function from Ben Bolker's page, linked above.

###########  CREATE LINK FUNCTION   ##################


logexp <- function(exposure = 1) {
  ## hack to help with visualization, post-prediction etc etc
  get_exposure <- function() {
    if (exists("..exposure", env=.GlobalEnv))
      return(get("..exposure", envir=.GlobalEnv))
    exposure
  }
  linkfun <- function(mu) qlogis(mu^(1/get_exposure()))
  ## FIXME: is there some trick we can play here to allow
  ##   evaluation in the context of the 'data' argument?
  linkinv <- function(eta) plogis(eta)^get_exposure()
  logit_mu_eta <- function(eta) {
    ifelse(abs(eta)>30,.Machine$double.eps,
           exp(eta)/(1+exp(eta))^2)
  }
  mu.eta <- function(eta) {
    get_exposure() * plogis(eta)^(get_exposure()-1) *
      logit_mu_eta(eta)
  }
  valideta <- function(eta) TRUE
  link <- paste("logexp(", deparse(substitute(exposure)), ")",
                sep="")
  structure(list(linkfun = linkfun, linkinv = linkinv,
                 mu.eta = mu.eta, valideta = valideta,
                 name = link),
            class = "link-glm")
}

###########



#####   Logistic Exposure   #####


###########
glm.dat%>%
  mutate(age=(day-exposure))->glm.dat # calculate age at the start of the exposure interval

mod.logexp<-glm(obs.ld~age,family=binomial(link=logexp(glm.dat$exposure)),
                data=glm.dat)

broom::tidy(mod.logexp)->parms.logit



parms.logit


#DSR
parms.logit%>%
  pull(estimate)->dsr.logit

age=seq(0,23,1)

dsr.logexp<-plogis(dsr.logit[1]+dsr.logit[2]*age) #Estimated daily survival rates

dsr.logexp  #Estimated DSR

S   # data-generating DSR values

#Comparing cumulative survival probabilities
prod(dsr.logexp)           
           
prod(S) 

And now how I think I am supposed to use the cloglog link + offset

glm.dat%>%
  mutate(age=(day-exposure),
         c.age=scale(age,scale=FALSE))->glm.dat # need to center age to avoid errors 


mean(glm.dat$c.age)->mu.cage


mod.clog<-glm(obs.ld~c.age+offset(log(exposure)),
              family=binomial(link='cloglog'),
              data=glm.dat)

broom::tidy(mod.clog)%>%pull(estimate)->parms.clog

age=seq(0,23,1)

c.age=age-mu.cage

#offset not included because log(1) (1 exposure day) = 0

lp= parms.clog[1]+parms.clog[2]*c.age # linear predicator

dsr.clog=1-exp(-exp(lp))

dsr.clog

And now to compare them all

#logistic exposure DSR
dsr.logexp

#cloglog DSR

dsr.cloglog

# "Truth"
S

#Cumulative survival probabilities
prod(dsr.logexp)
prod(dsr.clog)
prod(S)

I get why both methods might not match "Truth" due to sampling variability, but I was sort of expecting the logistic exposure and cloglog method to be much closer to one another.

Answer 1

The link functions can differ substantially:

I haven't checked your code (and coding questions per se are off-topic here), but if any probabilities of an event within a time period (given no event in a prior period) are reasonably high then you shouldn't be too surprised to find differences.

The cloglog link is appropriate under a proportional hazards assumption when events are observed only at the ends of discrete time intervals. It's sometimes called a "grouped proportional hazards" model.

Tutz and Schmid, in Chapters 3 and 4 of Modeling Discrete Time-to-Event Data, show some comparisons. If the underlying model truly follows proportional hazards with events only recorded at the ends of time intervals, "by definition the grouped proportional hazards model is correctly specified even if the interval lengths (and hence the number of ties) are large, RMSE values of this model are smaller than those of the Cox model in continuous time [with event times grouped]" (page 67). A continuation-ratio logit model performs about as poorly as the "Cox model in continuous time with event times grouped" if proportional hazards hold and time intervals are long.

Then again, if the data really followed a continuation-ratio logit model, a cloglog link would under-perform.

Code for plot:

curve(log(x/(1-x)),from=0.02,to=0.98,col="blue",bty="n",ylab="Link")
curve(log(-log(1-x)),from=0.02,to=0.98, add=TRUE,col="red")
legend("topleft","logit, blue\ncloglog, red",bty="n")