Seasonal survival using coxph

Question

I am having issues including "season" as a predictor variable in my survival analysis using cox proportional hazard models. I work with sage-grouse and have reason to believe their survival varies by season. I provided a sample data set below.

library(survival)
library(dplyr)

BirdID <- c(1,1,1,1,2,2,2,2,3,3,3)
Event <- c(0,0,0,0,0,0,0,1,0,0,1)
Start <- c(0,13,26,39,0,13,26,39,0,13,26)
Stop <- c(13,26,39,52,13,26,39,48,13,26,35)
Avg.Move <- c(2,5,6,9,15,9,16,25,1,3,18)
Season <- c("spring", "summer", "fall", "winter", "spring", "summer", 
"fall", "winter", "spring", "summer", "fall")

C <- cbind(BirdID, Start, Stop, Event, Avg.Move, Season)
T <- as.data.frame(C)
T$Start <- as.numeric(as.character(T$Start))
T$Stop <- as.numeric(as.character(T$Stop))
T$Event <- as.numeric(as.character(T$Event))
T$Avg.Move <- as.numeric(as.character(T$Avg.Move))

I have formatted my data so that each row is a different season because I have time-dependent variables (e.g. Avg.Move = average distance moved per season). I'm guessing that this may be part of my problem. Anyway, when I try to run basic coxph models with different combinations of my predictors (Season and Avg.Move) I continually get errors.

L <- coxph(Surv(Start, Stop, Event) ~ Season, data=T)
## Error in fitter(X, Y, strats, offset, init, control, weights = weights, : 
routine failed due to numeric overflow.This should never happen.  Please 
contact the author.

X <- coxph(Surv(Start, Stop, Event) ~ Avg.Move+Season, data=T)
## Warning messages:
1: In fitter(X, Y, strats, offset, init, control, weights = weights,  :
Loglik converged before variable  1,3,4 ; beta may be infinite. 
2: In coxph(Surv(Start, Stop, Event) ~ Avg.Move + Season, data = T) :
X matrix deemed to be singular; variable 2

I've tried coding Season as numeric (1,2,3,4) and binary (10000, 0100, 0010, 0001) but nothing has helped.

Thanks in advance,

Kyle

Answer 1

You only have three individuals and two events. Thus, you only have two partial likelihood terms on the log-likelihood scale. Thus, you are very limited by your data. In particular, you cannot estimate a dummy variable with four-levels as DJA mention's as you will have at-least two levels without any events. Thus, you can set the coefficient for these levels to minus infinity and get a hazard of zero (which is what you observe in the data but likely not what you expect).

The code below stress that you cannot estimate much from the data

da <- data.frame(
  BirdID   = c(1,1,1,1,2,2,2,2,3,3,3),
  Event    = c(0,0,0,0,0,0,0,1,0,0,1),
  Start    = c(0,13,26,39,0,13,26,39,0,13,26),
  Stop     = c(13,26,39,52,13,26,39,48,13,26,35),
  Avg.Move = c(2,5,6,9,15,9,16,25,1,3,18),
  Season   = c("spring", "summer", "fall", "winter", "spring", "summer", 
               "fall", "winter", "spring", "summer", "fall"))

# you only have three birds and two events
length(unique(da$BirdID))
#R [1] 3
sum(da$Event)
#R [1] 2

# you have hard time getting a good estimate of the non-parametric baseline
library(survival)
fit <- coxph(Surv(Start, Stop, Event) ~ 1 + cluster(BirdID), data = da)
plot(survfit(fit)) # plot survival curve

# You are still limited if you assume a constant intercept
fit <- survreg(Surv(Stop - Start, Event) ~ 1, dist = "exponential", data = da)
summary(fit)
#R 
#R Call:
#R survreg(formula = Surv(Stop - Start, Event) ~ 1, data = da, dist = "exponential")
#R             Value Std. Error    z       p
#R (Intercept) 4.212      0.707 5.96 2.6e-09
#R 
#R Scale fixed at 1 
#R 
#R Exponential distribution
#R Loglik(model)= -10.4   Loglik(intercept only)= -10.4
#R Number of Newton-Raphson Iterations: 5 
#R n= 11 

Answer 2

I believe the reason you are having an issue is because 2/4 levels of your factor "Season" have no events. Only winter and fall have a single event occur during those seasons. All other seasons haven 0 events occur. That is why you are not getting convergence.

See: Cliff AB (https://stats.stackexchange.com/users/76981/cliff-ab), Dealing with no events in one treatment group - survival analysis, URL (version: 2015-05-19): https://stats.stackexchange.com/q/153070