# Seasonal survival using coxph

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.

Kyle

• Hello, I have the same problem. did you ever figure this out? Jul 21 '17 at 16:57
• I guess you could try this: 1) table(T\$season)to is if it is singular; 2) sample size is only 11 obs, too small, try to get more obs and run it again. Sep 19 '18 at 9:42

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  3 sum(da$$Event)
#R  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


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